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■ In recent years the goal of estimating different cosmological parameters precisely has set new chal- 

lenges in the effort to accurately measure the angular power spectrum of CMB. This has required 
f"^ , removal of foreground contamination as well as detector noise bias with reliability and precision. Re- 

£SJ ' cently, a novel model-independent method for the estimation of CMB angular power spectrum solely 

from multi-frequency observations has been proposed and implemented on the first year WMAP data 
by Saha et al. 2006. All previous estimates of power spectrum of CMB are based upon foreground 
templates using data sets from different experiments. However our methodology demonstrates that 
CMB angular spectrum can be reliably estimated with precision from a self contained analysis of the 
WMAP data. In this work we provide a detailed description of this method. We also study and 
identify the biases present in our power spectrum estimate. We apply our methodoly to extract the 
power spectrum from the WMAP 1 year and 3 year data. 
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I. INTRODUCTION 
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Starting from the end of the last millenium remarkable progress in cosmology has been made by theprecise 
measurements of the anisotropies in CMB from different ground based as well as satellite observations P, H H, 0] • 
The Wilkinson Microwave Anisotropy Probe (WMAP) measures the CMB anisotropy over the 5 frequency bands at 
23 GHz (K), 33 GHz (Ka), 41 GHz (Q), 61 GHz (V) and 94 GHz (W). The observation system of the WMAP satellite 
^ ' consists of 10 differencing assemblies (DA), (ElEullII], one each for K and Ka bands, two for Q band, two for V band 
and four for W band. They are labeled as K, Ka, Ql, Q2, VI, V2, Wl, W2, W3, and W4 DA respectively. In the 1 
year and 3 year data release the WMAP science team has provided the science community with large amount of high 
quality data sets measured by these 10 DAs. However, extracting the primordial signal from these large data set is a 
non trivial task. The anisotropies in CMB are weak in comparison to those originating due to radiation in the local 
universe, which inevitably contaminate the observed signal. These dominant foreground emissions are from within 
the milky way as well as from extragalactic point sources In the low frequency microwave regime the strongest 
contamination comes from the synchrotron and free free emission [E E3] • At higher frequencies, where synchrotron and 
free free emissions are low, dust emission dominates. A reliable extraction of CMB signal from the multicomponent 
foreground contaminated data is thus complicated. There exists several methods in literature to remove foregrounds 
using foreground tracer templates [lol ITT1 . fl2| built from observations from other experiments. However, this requires 
; I ■ a prior model of spatio-spectral dependence for all the foreground components. The effect of uncertanities in the 
foreground models to estimate CMB anisotropies have been discussed in Refs. [ill 13 • The foreground cleaning is 
applicable in the region away from the galactic plane. Detector noise is another important concern that has to be 
addressed in order to precisely measure the angular power spectrum. The angular power spectrum from detector 
noise is dominant over the CMB power spectrum at the small angular scales. The detector noise, being a random 
quantity, is treated differently from foreground contamination which are treated here as fixed templates on the sky. 
However each of the DA's of WMAP has uncorrelated noise property [13, El, EE El- WMAP science team used this 
property to remove detector noise bias in a cross correlated power spectrum obtained from two different DA [ToL EH 
using DA maps with frequencies 41 GHz, 61 GHz and 94 GHz. 

An interesting model independent method to remove foregrounds from the multi-frequency observations of CMB 
without any assumption about foreground components has been proposed in Ref. [l7| and implemented in Ref. [Hj] 
in order to extract the CMB signal from the WMAP data. The foreground emissions were removed by exploiting 
the fact that their contributions in different spectral bands are considerably different while the CMB anisotropy 
power spectrum is same in all the bands in unit of thermodynamic temperature due to the Planck blackbody energy 
spectrum of the CMB, [TEHS]. The main advantage of this foreground cleaning method is that it is totally free from 
any assumption about foreground modeling. Another advantage is that it is computationally fast. However the auto 
power spectrum obtained from a single cleaned map as reported in Ref. [1 81 ] is not directly usable for primordial power 
spectrum estimation at smaller angular scale. This is because the detector noise bias dominates over the CMB power 
spectrum at smaller angular scale, beyond the beam width of WMAP detectors. 

In earlier publications [U, [22|, H3] we extended this model independent foreground removal method to remove 
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detector noise bias also. In this work we describe the basic formalism of our previous work in detail. We apply our 
method both on the WMAP 1 year and WMAP 3 year data to estimate CMB angular power spectrum. We form 
several cleaned maps using different cross combinations of the DA maps. Finally we form several cross power spectra 
where detector noise bias is removed. The results from this analysis are summarized in figure In this figure we 
show the power spectra estimated from the WMAP's 1 year and 3 year data using the multi- frequency combination 
of CMB maps. The bottom part of this figure shows the residual unresolved point source contamination that were 
corrected for to obtain these two spectra. These power spectra are obtained without making any explicit model of the 
foregrounds or the detector noise. In most power spectrum extraction procedures, only the three highest frequency 
channels observed by WMAP have been used to extract CMB power spectrum. We present a more general procedure 
where we use observations from all the five frequency channels of WMAP. The primary merit of the foreground removal 
method is that it avoids any need to remove foregrounds based upon extrapolated flux measurements at frequencies 
far away from observational frequencies of WMAP [l(J HH HH . 

Presence of a bias in the internally cleaned maps has been reported earlier in the Refs. [HI, [2(| HtJ . In this work we 
also perform a detailed study of the nature of bias in the cleaned power spectrum. We show that the cleaned power 
spectrum is not an unbiased estimator of the underlying CMB spectrum. Naively one expects that there might be 
some residual foreground contamination causing a positive bias in the power spectrum. For a simplified approach, 
which cleans the entire sky simultaneously, without sub-dividing into regions of varying foreground contamination, 
we are able to analytically compute the cleaned power spectrum in terms of a CMB signal and the foreground plus 
detector noise covariance matrix. The existence of bias is easily identified from these results. We also report and 
quantify an interesting negative bias in the cleaned power spectrum. This negative bias is directly determined by 
the underlying CMB power spectrum and is strongest at the lowest multipoles. The analytical results for the bias 
estimations make the model independent foreground removal method more interesting for use in cosmology. The 
bias in the cleaned power spectrum can be removed following models of foreground and noise covariance matrices 
after a model independent foreground removal is performed. This makes the estimated power spectrum less prone to 
uncertainties of the foreground modeling compared to a method which tries to minimize foreground using external 
templates. In this work we debias the CMB anisotropy spectra obtained from WMAP data only at the large multipolc 
regime, I > 400. It turns out that because of large noise level of WMAP the point source bias is the only issue in this 
multipole range. We do not attempt to employ a debiasing method at the low multipoles where the negative bias is 
expected to dominate, because of its complicated nature in the case of the iterative, multiregion cleaning method. A 
more detailed study on this issue will be reported in a future publication. 

The plan of this paper is as follows. The basic formalism and methodology to obtain power spectrum is described 
in the section [TTJ Here we also discuss the bias present in this method. The implementation of our methodology 
on the WMAP I year and 3 year data is discussed in section HTT1 Here wc also obtain an analytic expression for the 
residual unresolved point source power contamination. The results are described in section ITVl and finally we conclude 
in section [V] Some of the notations used in this report are as follows. We denote matrices and vectors by boldfaced 
letters. Any variable (scalar, vector or matrix) which depends on the stochastic component is denoted by a hat on 
the top. As an example we note that the CMB power spectrum is a stochastic variable which is represented as Cf. 

II. BASIC FORMALISM AND METHODOLOGY 

In this section we outline the mathematical formalism to obtain the underlying power spectrum. We also quantify 
the bias present in the estimated power spectrum using our procedure. 

A. Foreground Removal 

The observed signal at frequency channel i in a differential telescope like WMAP can be modeled as 

AT l {h) = J (AT c (rV) + AT u (h')) B\h ■ h')dh' + AT m (n) . (I) 

Here AT c (h) and AT H (n) are respectively CMB and foreground component of the anisotropy in the channel, i. The 
detector noise in the channel, i, is AT n \ The beam function B l (h ■ n r ) represents the smoothing of the map due to 
finite resolution of the antenna of channel, i. The beam is assumed to be circularly symmetric as done in most analysis. 
We note that the detector noise is not affected by the beam function. An experiment such as WMAP provides multi- 
frequency maps AT l (h), i — 1,2,..., n c corresponding to observation of the CMB at n c different frequency bands. 
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Equivalently, in the spherical harmonic representation 



ai m = la 



lm 



Hra ■ 



(2) 



where a; m are respective spherical harmonic contributions to the maps and B\ are Legendre transform coefficients of 
the beam B l (n ■ n'). The aim is to linearly combine the maps with appropriate weights to get an optimal estimator 
of the CMB anisotropy AT c (n') that minimizes the contribution from AT f (n'). The linear combination of the multi- 
frequency maps available can be carried out in pixel space, or, in the equivalent representation in terms of the spherical 
harmonic coefficients. The former has been followed by the WMAP team to produce Internal Linear Combination 
(ILC) map and also in a related, but more elaborate approach in Ref. [27] to produce LILC map. The approach of 
carrying out a multi- frequency maps combination in the spherical harmonic space was proposed in Ref. [17]. For the 
first year of WMAP data this was implemented in Ref. [IH . This method has the advantage that we can simultaneously 
take into account variation of foreground with sky positions and with different multipoles for a given sky position. 
We define a cleaned map as a linearly weighted sum of the maps at different frequencies, 
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Here W[ is a weight factor which depends upon the multipole I and the frequency channel, i. Since each of the 
channels has a different beam resolution, the maps are deconvolved by the corresponding circularly symmetric beam 
transform functions Bi prior to the linear combination. The total power in the cleaned map at a given multipole I is 
then 
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Substituting eq. ([3]) in eq. Q we obtain, 
where the matrix Ci is given by 



(jClean = W^Wf 
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Wi is a row vector describing the weights for different channels, 

W 1 = (wl 7 wf 



1 ,uj 1 ,,...,W 1 
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and C l is the cross power spectrum between the i and j channel, 
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By construction, the Ci matrix is symmetric. The spherical harmonic coefficients extracted from a map contain CMB 
signal and foregrounds, both smoothed by the beam function of the optical instrument used in the experiment, as 
well as detector noise. Using eq. ^ we obtain, 
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where 
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is used to denote the total non-CMB contamination in the map. Since the CMB contribution is indepen- 



dent of frequency, the expression for B i B j simplifies to 
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for all values of i and j. Now using eq. ([5]) we obtain 

fjClean = (jS^^T^T + WjCp"^ , (10) 

where eo is a column vector with unit elements 

eo = | ;; | . (ii) 

The CMB signal power in the cleaned map is kept unaltered by imposing the constraint 

W ieo = e^W T = 1 , (12) 
on weights W\. Using eqs. (flTJ|) and (fl2|) . the total power at a multipole I in the cleaned map is 

£Clean = £S + W^^W? . (13) 

Thus the CMB signal power is only an additive positive constant in the expression of the total power of the weighted 
map. Hence, choosing weights that minimize (7 ; clcan also minimizes the combined contamination coming from fore- 
ground and detector noise. It may be shown easily, using the Lagrange's multiplier method, that minimizing Cp lean 
subject to the constraint eq. p"2|) . gives the following expression for the weight factors [l7l |26|. 

e T C 1 

w. = -S3r- • ( 14 ) 

eg C, e 

Following eqs. (0) and (fT4)) we can express the power in the cleaned maps neatly as 

£Clean = * . (15) 

e C, e 

It is important to note that in cases when Ci is singular, it is possible to generalize the above expressions for weights 
and cleaned power spectrum in terms of the Moore-Penrose Generalized inversion (MPGI) of Ci. The generalized 
weights and cleaned power spectrum are then given by 

W. = , (16) 

Cf lean = — if- . (17) 
For further details of the analytic derivation of the weights and cleaned power spectrum we refer to appendix [SJ 

B. Biases in the foreground cleaning method 

Although the foreground cleaning is performed with the constraint that CMB power spectrum remains preserved, 
the method biases the final power spectrum. This section is devoted to a discussion of the full bias in the method. 

The existence of some bias is not difficult to anticipate and understand. The method is intended to perform a 
minimization of the foreground power spectrum which is a positive definite quantity. Unless the foreground cleaning 
is fully effective at all multipoles, the minimization would leave some residual foreground which naively would give 
rise to a positive bias in the cleaned power spectrum. However, it is very interesting that there exists an additional 
negative bias in the method. This negative bias is strongest at the lowest multipoles and increases in magnitude with 
increase in number of channels that are combined. 

Let us consider n c number of channels in the linear combinations for the cleaned map. The maps at each channel 
consists of the CMB signal and foreground contamination coming from the different foreground components (syn- 
chrotron, free- free, dust etc). Although, the number of distinct foreground components could be arbitrary, for brevity 
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and simplicity, here we club the contributions from all the components into a single foreground term. The foreground 
contribution can be described by a covariance matrix. Later, after simplification of the expressions, we can split up 
the total foreground covariance matrix in terms of the covariance matrices of the distinct constituent components. In 
the entire discussion that follows we do not treat foregrounds as a stochastic component. Instead we consider them as 
fixed templates in all the realizations. This is entirely justified. We are simply interested in computing a foreground 
free template rather than estimating information about the distribution of foreground from which they are drawn. 
We discuss the bias in few different cases depending upon the rank of the covariance matrix Ci. 

1. Case: Rank Ci < n c — 1 

First we consider an ideal case wherein detector noise remains absent and each foreground component follows a 
rigid frequency scaling on the entire sky. Let the number of foreground components be n/. Then the rank of the 
foreground covariance matrix, Cf, is also n/. 

The p th foreground component for channel i is denoted by Fq(8, 4>)fp, where the frequency dependence, /* and the 
spatial (sky) dependence F p (8, </>) are explicitly separable in the rigid scaling assumption. (Here, Fq(9,(/)) is the p th 
foreground template based on frequency i^o, so that = 1, for frequency i/q.) We denote the CMB component by 
C(8, 4>). Full signal map at frequency channel, i, is then given by 



S\QA)=C(9A)+Y. F ^ e ^)f l v 



(18) 



Alternatively, in the spherical harmonic space, 
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The auto power spectrum of the i channel 
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In the above equation, C, ^° is the correlation between any two foreground components p,p' and C^ JU> '" denotes 
the chance correlation between CMB signal and p th foreground component. The cross power spectrum between two 
channels i, j is given by 



c/(p)0 



C\ J = Cf 



It is convenient to define the vectors, 
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(22) 



fi° = £f> p0 - 

P =i 

A little algebraic manipulation allows us to recast eq. pTjl as 

C, = CfeoeJ + f,°ej + e t° T 



(23) 



(24) 



The above equation will be useful to compute the bias in the cleaned power spectrum. On the ensemble average the 
cleaned power spectrum as given by eq. (|17p could be simplified with the help of successive use of a set of theorems 
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reported in Refs. [23, l23| • An elaborate discussion of these theorems is given in appendix [E] Assuming statistically 
isotropic CMB sky we obtain the following expression for the ensemble average of the cleaned power spectrum 

(Cf^) = (df)-n f ^ l . (25) 

We easily note a few interesting aspects of the above equation. First of all there exists a negative bias. Secondly, 
because of ~ 1/(2Z+ 1) decay, this bias is important at the lowest multipoles. Another important point to note is that 
the bias depends on the underlying CMB power spectrum. Therefore it is possible to debias any given statistically 
isotropic CMB model in this approach by constructing an appropriately scaled estimator. Lastly, we find that there is 
no foreground bias. One would naively expect the foregrounds would contribute a positive bias in the power spectrum. 
However this does not happen in this case as the rigid scaling assumption along with the condition n c > rif + 1 ensure 
that sufficient amount of spectral information are available to remove all the foregrounds. The negative bias arises 
as the weights are to be determined from the empirical covariance matrix to take into account information available 
from the observed data. 



2, Case: Rank (Cf) = n c 

The rigid scaling assumption for the foreground contaminants considered in the previous section is at best a 
reasonable approximation and is known not to be valid in general. As mentioned in Ref. @ a foreground component 
with varying spectral index over the sky could be approximated in terms of two templates, provided the variation 
is small compared to the mean spectral index over the sky. A stronger variation will need more than two templates 
for reasonable modeling. In such a situation, if the number of templates required for modeling of all the foreground 
components exceeds the number of maps available for linear combination then Cf is of full rank. In this case a positive 
foreground bias appears along with a negative bias. The negative bias is similar to the previous case in that it remains 
(c a ) 

proportional to ■ Keeping aside the detector noise for the moment, the ensemble average of the cleaned power 
spectrum is given by 

(Cf lean ) = (6f) + L_^ + (i_ nc )l5_) (26) 

\ / \ V e T(C[)" 1 e ' '21 + 1 

Detailed derivation of the above equation is similar to derivation of eq. ([25]) . As Cf is of full rank and a positive 
definite matrix, the second term on the right causes a positive foreground bias. 



3. Noise and foreground case 

The discussion in the last two sections does not consider any detector noise. We now consider the most general case 
where we have both foreground and detector noise. Following a method similar to that used in derivation of eqs. (j2"5|) 
and (|26[) we can show that, on the ensemble average the cleaned power spectrum in this case is given by 

/ £Clean\ =/&) + / — 

V ' V ' W(Cf+ N ) 

We have carried out Monte-Carlo simulations to verify the analytical result given by eq. (|25|) . We perform simulations 
of foreground cleaning using n c = 3 channels, corresponding to 41,61 and 94 GHz frequencies. The foreground 
model consists of a synchrotron and free-free emission. Each of the foreground component is assumed to follow a 
rigid frequency scaling all over the sky. First we generate synchrotron and free-free templates using the Planck Sky 
Model 1 at the 3 different frequencies. Since these templates do not follow rigid frequency scaling all over the sky, ( 



+ (1 - n c ) 



21 + 1 



(27) 



1 We acknowledge the use of version 1.1 of the Planck reference sky model, prepared by the members of Working Group 2 and available 
at www.planck.fr/heading79.html. 
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eg., synchrotron spectral index varies with the position of the sky), we developed a method to regenerate rigid scaling 
synchrotron and free templates at these 3 different frequency channels 2 






Input C C | - Biased Cleaned C, 




2C C |/(2I+1) 







100 200 
Multipole, I 



100 200 
Multipole, I 



FIG. 1: The negative bias in the extracted power spectrum at low I is shown by the red line in the left panel. The bias corrected 
spectrum is plotted in blue line. However it lies entirely behind the green line and is not visible. The right panel shows how 
well the analytic results for bias match with those obtained from Monte Carlo simulations. 



The left hand panel of the fig. [T] based on 1000 Monte Carlo simulations shows that there exists a negative bias 
in the cleaning method. For the assumed model of foreground components with rigid frequency scaling, the second 

term on the right in eq. (|25[) contributes —2^Cf^/2l + 1 as a negative bias. In the right hand panel of the fig. Q] 

(c c ) 

we explicitly show that the magnitude of the negative bias is exactly compensated by 2 foH ■ The negative bias is 
important at the lowest multipoles and becomes negligible at high I, e.g. the bias is only w 3/j,K 2 at multipole I 
« 800. The bias corrected spectrum plotted in green color in the left panel of this fig. Q]is completely hidden by the 
blue curve corresponding to the input CMB power spectrum which is used to generate random realizations of CMB 
maps. 

The bias in the cleaned power spectrum for different cases described so far is computed following the assumption 
that foreground cleaning is simultaneously carried out over the entire sky. However for the cleaning purpose of 
WMAP maps we follow a more sophisticated method. The main reason behind opting for a sophistication is that the 
spectrum of foreground emission as well the amplitude strongly depend on the location in the sky. The effectiveness 
of the foreground cleaning in varying foreground spectral index has been studied in the existing literature, 0, 
Following the discussion in section III B 21 a varying foreground spectra would cause a positive foreground bias in the 
cleaned power spectrum. To minimize this bias we followed Ref. [l8T ] and divide the sky into 9 regions based on a 
simple estimation of the level of foreground contamination (see appendix |5]). We then carry out foreground cleaning 
iteratively taking one region at a time starting with dirtiest region. The advantage of such iterative method is that 
diffuse foreground contamination which is dominant at the low multipoles compared to the detector noise could be very 
effectively minimized without a precise model of residual bias. The main principle for iterative cleaning is to search 
for sky regions where foreground emission could be assumed approximately constant or at least nearly constant. Such 
a philosophy of foreground cleaning has been adopted by WMAP team to produce their Internal Linear Combination 
Map In a future publication we would generalize the bias results reported in this article to the most general 

method of iterative cleaning. Another advantage of iterative method is that one has the freedom to adjust weights on 
different sky positions depending upon the amplitude of the foreground components leading to a better cleaning. 



2 For this purpose, we first note that for each component the approximate normalization yC ; l /C ; 23 remains roughly constant for all I. 
Here C 23 is the synchrotron (or free -free) power spectrum for frequency 23 GHz and C ; 8 is the synchrotron ( or free-free) power spectrum 
at any of other 3 frequencies. We generate a synchrotron ( or free-free ) template at the i th frequency channel following a scaling of 
the 23 GHz template by the number ^Cf/Cf 3 . This ensures that the emisssion for each foreground component at different frequencies 
follows rigid frequency scaling. 
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n c -channel combination 


Cleaned Map 


Dj+Df + . 

+ D| + . 
DJ + D| + . 


. + D" C 
. + D™<= 
. + D" C 


Ci 

c 2 
c 3 


+ Dl + . 


, + D- 


c d 



TABLE I: List of several possible cleaned maps with uncorrelated detector noise properties. Here we assume n c number of 
frequency bands with each band having d number of independent detectors. 

C. Power spectrum Estimation 

The power spectrum is estimated by cross correlating multichannel foreground cleaned maps from independent 
Differencing Assemblies (DA) . The alternative is to use auto correlation and use a detector noise model to debias the 
power spectrum. In this case the detector noise model uncertainty affects the mean of the cleaned power spectrum. 
We note that the WMAP team too used cross-correlation of independent DA to obtain their power spectrum. Let 
us consider an hypothetical experiment with n c number of frequency channels each of which is denoted by i. Each 
frequency channel consists of d number of independent detectors denoted by D*. Thus D* is the map observed by 
j th detector of i th channel. We then propose to form several cleaned maps in such a way that each cleaned map 
consists of one map from each of the n c number of channels. By choosing appropriate combinations of DA maps we 
identify two different cleaned maps having totally disjoint set of detectors. Then assuming that the noise properties 
are uncorrelated in two different detectors, the noise bias will not affect the cross power spectrum of two such cleaned 
maps at the I s * order. The final power spectrum could be a simple average of all possible such cross power spectra. 
In table Q] we show a set of possible cleaned maps Cj,j — 1, 2, .., d with uncorrelated detector noise properties. Hence 
number of cross power spectra can be obtained that do not have any noise bias. In general the number of 
possible cleaned maps is more than d. If we divide the combination of cleaned map 1 in p number of subsets and the 
combination of cleaned map 2 in q number of subsets then one only needs to make sure that these p and q subsets 
have uncorrelated noise properties. 

III. IMPLEMENTATION ON WMAP DATA 

In our method we linearly combine maps corresponding to a set of 4 DAs at different frequencies. 3 We treat the K 
and Ka maps effectively as observations of the CMB sky in two different DA at the low frequencies. Therefore, we use 
either K or Ka maps in any combination. In case of the W band, 4 DA maps are available. We simply form pairwise 
averaged map taking two of them at a time and form 6 effective DA maps corresponding to W band. Wij represents 
simply an averaged map obtained from the i th and j th DA of W band. In table |TT] we list all the 48 possible linear 
combinations of the DA maps that lead to 'cleaned' maps, Ci and CAi's, where i = 1, 2, . . . , 24. In an alternative 
approach, we form combinations to form cleaned maps excluding the most foreground contaminated K and Ka bands 
(referred to as the three-channel n c = 3 case). This leads to a total of 24 cleaned maps. All such combinations are 
shown in the bottom panel of the table HT1 In this combination the cleaned maps are labeled as Ci s, where i = 1, 2, 
•■■,24. 

The entire method leading to the power spectrum estimation consists of three main steps. First we perform 
foreground cleaning using several multi-channel combinations of WMAP maps. At the second step we obtain cross 
power spectra from these foreground cleaned maps. The foreground cleaning is similar to Refs. Finally we 

correct for the estimated residual unresolved point source contamination. We note in passing that each of these steps 
are logically modular and each of them could be modified or improved independent of the other. 



3 In Ref. [lS( . a single foreground cleaned map is obtained by linearly combining 5 maps corresponding to one each for the different 
WMAP frequency channels. For the Q, V and W frequency channels, where more than one maps were available, an averaged map was 
used. However, averaging over the DA maps in a given frequency channel precludes any possibility of removing detector noise bias using 
cross correlation. 
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4-channcl combinations (n c — 4) 




(K,KA)+Q1+V1+W12=(C1,CA1) 
(K,KA)+Q1+V1+W13=(C2,CA2) 
(K,KA)+Q1+V1+W14=(C3,CA3) 
(K,KA)+Q1+V1+W23=(C4,CA4) 
(K,KA)+Q1+V1+W24=(C5,CA5) 
(K,KA)+Q1+V1+W34=(C6,CA6) 
(K,KA)+Q2+V2+W12=(C7,CA7) 
(K,KA)+Q2+V2+W13=(C8,CA8) 
(K,KA)+Q2+V2+W14=(C9,CA9) 
(K,KA)+Q2+V2+W23=(C10,CA10) 
(K,KA)+Q2+V2+W24=(C11,CA11) 
(K,KA)+Q2+V2+W34=(C12,CA12) 


(K,KA)+Q1+V2+W12=(C13,CA13) 
(K,KA)+Q1+V2+W13=(C14,CA14) 
(K,KA)+Q1+V2+W14=(C15,CA15) 
(K,KA)+Q1+V2+W23=(C16,CA16) 
(K,KA)+Q1+V2+W24=(C17,CA17) 
(K,KA)+Q1+V2+W34=(C18,CA18) 
(K,KA)+Q2+V1+W12=(C19,CA19) 
(K,KA)+Q2+V1+W13=(C20,CA20) 
(K,KA)+Q2+V1+W14=(C21,CA21) 
(K,KA)+Q2+V1+W23=(C22,CA22) 
(K,KA)+Q2+V1+W24=(C23,CA23) 
(K,KA)+Q2+V1+W34=(C24,CA24) 


3-channcl combinations (n c — 3) 




Q1+V1+W12=C1 
Q1+V1+W13=C2 
Q1+V1+W14=C3 
Q1+V1+W23=C4 
Q1+V1+W24=C5 
Q1+V1+W34=C6 
Q2+V2+W12=C7 
Q2+V2+W13=C8 
Q2+V2+W14=C9 
Q2+V2+W23=C10 
Q2+V2+W24=C11 
Q2+V2+W34=C12 


Q1+V2+W12=C13 
Q1+V2+W13=C14 
Q1+V2+W14=C15 
Q1+V2+W23=C16 
Q1+V2+W24=C17 
Q1+V2+W34=C18 
Q2+V1+W12=C19 
Q2+V1+W13=C20 
Q2+V1+W14=C21 
Q2+V1+W23=C22 
Q2+V1+W24=C23 
Q2+V1+W34=C24 



TABLE II: The table on the top shows 48 different combinations of the DA maps used in our 4 channel cleaning method. The 
bottom table shows the 24 possible combinations in the 3 channel cleaning method. 



A. Map cleaning 

We use the 'raw' DA maps (i.e., that have not undergone any foreground cleaning process) both for the WMAP 1 
year and WMAP 3 year data release from the LAMBDA website. These maps follow HEALPix 4 pixelization scheme 
at a resolution level N s i^ e = 512 corresponding to approximately 3 million sky pixels. All these maps are provided 
in the 'nested' pixelization scheme which is suitable for nearest neighbor searches. However for converting maps 
to spherical harmonic space and vice versa it is computationally advantageous to convert them to ring format that 
facilitates the use of Fast Fourier transformation method along equal latitudes. Therefore prior to the analysis all the 
maps corresponding to 10 DA's were converted to 'ring' pixelization scheme. When converting a map to spherical 
harmonic space we restricted ourselves to a maximum multipole, l max = 1024. 

The spectrum of foreground emission has some dependency on the location on the sky. A better cleaning may 
be achieved if we partition the entire sky into certain number of sky parts depending upon the level of foreground 
contaminations [1 81 ] . Then cleaning is done for each sky-parts iteratively. For each regions we will have then different 
weights which are chosen to minimize foreground contamination from that particular region. Cleaning becomes more 
efficient by allowing the weights to depend on the sky parts rather than using a single set of weights for the entire 
sky. Following Ref. [Tc| we partition the sky depending upon the level of foreground contamination. 5 

The entire sky is partitioned in a total of r = 9 regions depending upon their foreground contamination. All the r 
sky parts are shown in the fig. [2j Individual sky parts are color-coded according to the number r assigned to them. 
The dirtiest part is labeled with the maximum index. The procedure of forming these sky parts are similar to Ref. [LSj 



4 For comprehensive studies about HEALPix we refer to Ref. 30, 31, 32]. 

5 Although this scheme is found to be effective, it is possible to envisage other schemes such as those that make foreground contamination in 
each part closer to rigid scaling approximation. In an ongoing study we have investigated and demonstrated improvement in foreground 
cleaning for different levels of partitioning [33ll . 
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FIG. 2: The 9 different masks in orthographic projection that correspond to the partitions of the sky. Each region is coior 
coded to distinguish them visuaily. The cieanest region is the region farthest away from the galactic plane. The dirtiest region 
(black spots) lie in the galactic plane. 



with small variants. We describe our sky partitioning to identify these r regions in detail in appendix IB"1 

For each combination, the cleaning procedure is performed in r iterations starting from the dirtiest region. As 
in Ref. [l8j we call the initial foreground contaminated maps as the initial temporary maps. At the end of each of 
iteration we obtain a set of partially cleaned temporary maps which are used as input in the next iteration. In i th 
iteration (i = 1, 2, 3, r) we perform the following three steps : 

1. Using the i th mask we select the i th region of the sky from the temporary maps. Then we obtain power spectrum 
matrix from the i th region only. 6 

2. Using the weight factors we obtain a full sky cleaned map using eq. ([3]). 

3. Next we replace the z th region of the temporary maps by the corresponding region of the cleaned map. Before 
replacement the cleaned map is smoothed to the resolution of the different frequency channels. After replacement 
the maps define new temporary maps for the next iteration. 

The entire cleaning procedure is automated to carry the three steps iterated r times to obtain the final cleaned 
maps. The final cleaned maps have the resolution of the W band. 

In eq. |7]) the weights could have numerical errors if the Ci matrix becomes ill-conditioned at specific values of 
I. This can happen if by chance any mode has almost equal contribution combined from the CMB, foreground and 
detector noise. Hence, in practice, the numerical implementation obtains weights using the Ci matrix smoothed over 
a range of Al = 11 prior to inversion. 

As shown in table [Til f° r the four channel combination (n c — 4), a total of 48 cleaned maps can be obtained using 
all the possible combinations of the DA maps for each of the WMAP 1 year and WMAP 3 year data. The cleaned 
map C8 for WMAP 1 year data is shown in the left panel of the figure [3J There is some residual foreground left in 
the galactic plane as seen in this map. We apply Kp2 mask, supplied by the WMAP team, to flag the contaminated 
pixels near the galactic plane. Therefore these residuals do not affect our final estimated power spectrum. A similar 
map is obtained from the WMAP 3 year data also. The right panel of the figure [3] shows CA8 map from WMAP-3 
data. Both the maps look similar. 

The weights Wi are shown in the figures Q] for one of the cleaned map (CA1) for the cleanest and second cleanest 
region. For low I where the diffuse foregrounds are dominant, the weights take large positive and negative values to 



These power spectrum matrix is obtained from the partial sky spherical harmonic coefficients. In practice while using 1 year WMAP 
data we have not combined K band map for I > 603 (Bi < for K band for I > 603 and the negative values of beam function specified is 
unphysical.). KA and Q band beam functions are available till I = 850 and / = 1000. Therefore we do not combine these bands beyond 
these values. Similar limits determined by the condition Bi > were used for 3 year data analysis also. 
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FIG. 3: The left panel shows the cleaned map C8 using the 1 year WMAP data. There is some residual foreground contami- 
nation near the galactic regions. The right panel shows the CA8 map using the 3 year WMAP data. The temperature scale of 
both the figures are chosen from the sky part outside the Kp2 mask. 
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FIG. 4: The left panel of the figure shows the weights for the cleanest region for the combination KA, Q1,V1,W12. The right 
panel show the weights for the same combination but for the second cleanest region. At low multipole foregrounds dominate. 
Therefore, weights take large positive and negative values to subtract foregrounds. 



subtract the foregrounds. At large I the maximum weight is given to the W band channel since it has the highest 
resolution. 



B. Power spectrum estimation 

Our final power spectrum is based upon the MASTER estimate introduced in Ref. [34]. Even after performing 
foreground cleaning the galactic disk region remains significantly contaminated. Hence one needs to exclude this 
region before estimating the cosmological power spectrum. However, flagging the contaminated sky region effectively 
introduces uneven weighting of the pixels. Those flagged have efectively zero weight and rest have unit weight. 
Weighting a map in the pixel space causes neighboring multipoles to get correlated. The correlation of the multipoles 
are described by the coupling matrix which is determined by the geometry of the sky coverage. Moreover, the power 
spectrum of such a map is biased low because of less sky coverage. The MASTER method debiases the partial sky 
power spectrum on the ensemble average by inverting the coupling matrix. 

The MASTER method originally implemented in Ref. [34| deals with auto power spectrum only. This can readily be 
generalized to the case of cross power spectrum [35[ . We first exclude the galactic region from all the 48 cleaned maps 
(Table |Hln c = 4) using Kp2 mask as supplied by the WMAP team. Then map pairs Ci & Cj to be cross correlated 
are chosen such that they do not have any DA common between them. This choice ensures that the detector noise 
bias does not affect any of the cross power spectrum. A total of 24 cross power spectra could be obtained which 
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are unaffected by the detector noise bias. Each of these cross spectra are then debiased from the partial sky effect 
following Ref. [35} using the coupling (bias) matrix corresponding to the Kp2 mask. The small scale systematic effects 
of beam and pixel smoothing were removed using appropriate circularized beam transform (34| and pixel window 
functions as supplied by the HEALPix package. The 24 cross power spectra are then combined with equal weights 
into a single 'Uniform average' power spectrum. 7 The final power spectrum is binned in the same manner as the 
WMAP's published result for ease of comparison. 

Both one and three year power spectra were obtained following this method. We defer a more detailed discussion 
on the 3 year power spectrum results to the section IfVBI We also estimate a residual contamination in the 'Uniform 
average' power spectrum for both WMAP 1 year and WMAP 3 year data from the unresolved point sources. A point 
source model used by WMAP team [l(| E2] an d estimated entirely from the WMAP data [ll[ is sufficient for our 
estimation. We recover a flat, (approximately 140/j.K 2 for WMAP 1 year data and 100/iK 2 for WMAP 3 year data), 
residual point source contamination in the two 'Uniform average' power spectra for I > 400. This residual is much less 
than actual point source contamination in Q, KA or K band and intermediate between V and W band point source 
contamination. 



C. Estimation of residual unresolved point source power spectrum 

There have been several studies regarding the point source contamination in the CMB maps [lH, H3, [HI]. We 
estimate residual contamination due to unresolved point sources in the 'Uniform average' power spectrum following 
the point source model constructed by the WMAP team [l(| Ojj]. The WMAP point source model consists of a 

point source covariance matrix in thermodynamic temperature unit following Cf s ^ — cuj)A (vi/vo) 2 { v j/ v o) 2 f° r 
two different frequency channel i and j. Here A is the amplitude of unresolved point source spectrum in antenna 
temperature at a reference frequency v$ and Cy is the conversion factor from antenna to thermodynamic temperature. 
We compute the residual unresolved point source contamination in each of the 24 cross power spectra assuming 
that the residual (extra galactic) point source contamination has a statistically isotropic distribution in the sky. As 
described in section Ull A[ the foreground cleaning procedure consists of r number of iterative steps. The total point 
source contamination in any of the cross power spectra consists of contribution from each of the individual sky parts. 
The final residual unresolved point source contamination is the uniform average over all such 24 residual point source 
contamination. Below we describe in detail residual point source contamination in a cross combination of two maps 

As stated in section Ull Bl we apply Kp2 mask on the cleaned maps prior to cross power spectrum estimation. Apart 
from removing the galactic region, Kp2 mask also removes a circular region of 0.6 degree radius around each of the 
resolved point sources. Therefore the resolved point sources cannot affect the cross power spectra. However in the 
detector noise dominated region, I > 400, the iterative foreground cleaning method cannot remove all power from 
the unresolved point sources. As resolved point sources are already masked it is important to estimate only residual 
unresolved point source contamination in each of the cross power spectrum. 

Let us assume that a\ m is the spherical harmonic coefficients obtained from cleaned map Ci, where i = 1, 2, 48, 

a\ m =B\at m + B\a^ . (28) 

Here a%}^ is the residual unresolved point source contamination. We note that we have not considered any detector 
noise contribution in this equation. However this does not imply any loss of generality in the discussion of cross 
power spectrum. The noise bias does not affect the cross power spectrum. Also we do not include any residual 
diffuse foreground in our study in this section. At the large multipole region where point source contribution becomes 
significant diffuse galactic contaminations are entirely subdominant and thus they are not a concern. 

We assume a statistically isotropic model of the residual unresolved point source distribution over the sky. Point 
sources are uncorrelated with CMB. On the ensemble average a partial cross power spectrum obtained from cleaned 
maps i and j is related to the CMB and point source power spectra as follows 

(Cf) = M W ((Ct) + (Cf»)) B\,4pf, . (29) 



7 There exists the additional freedom to choose optimal weights for combining the 24 cross-power spectra which we do not discuss in this 
work. 
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Here, Mw is the coupling matrix, Cf,^ is the estimate for the full sky residual unresolved point source power spectrum 
and B\,B\,pf, denote combined effect of beam and pixel smoothing. We can recast the above equation as, 



^-^)- <30) 

The next task is to obtain the estimates Cf, themselves in each cross spectra. To compute them, let us assume that 
f l (9, (f>) is the residual point source function present in the \ th cleaned map. Our cleaning method partitions the sky 
into 9 parts. The entire sky is then cleaned in a total of 9 iterations. If gk(9,4>) is the point source residual present 
in the k th sky part, we have 

k=9 

f(M)=E^(^)- (31) 

k=l 

After expanding both sides in spherical harmonics we obtain 

fc=9 

«L = E^»- ( 32 ) 

k=l 

The partial sky unresolved residual point source modes could be written in terms of the full-sky modes using 

®lm = E ^hnl'm' a l'm' • (33) 
I'm' 

We note that, af, m , represents the residual point source contamination in the entire cleaned map obtained after the 
k th iteration. This is different from the contamination, 5; TO , that is actually present in the k th sky part. The symbol 
Ml Li/ m / denotes the mode mode coupling matrix for the spherical harmonic modes for the k th sky part. We can 
rewrite af, , in terms of the spherical harmonic coefficients of the temporary maps obtained at k th iteration as, 

»=4 

4m'=J2 a l^V, (34) 
i=l 

where wf l is the weight for the i th channel at the k th iteration for multipole I. Following eqs. (|3^|) . (f3"3")l and (fM)) we 
obtain 

k=9 i=4 

4m = E E ^hnl'm' E a V m '^ ■ (35) 
k=l I'm' i—1 

The residual unresolved point source spectrum in one cross combination (i, j) is given by 

rn — l I i *j i 
^Ps(ij)\ _ \ a lm a lml 

^ 21 + 1 



\ c i )= E ^^rvf 1 - ( 36 ) 

Substituting aj m and a,^ we find, 

m—l k 7 k'—9 i,i'=4 

(^ S(ij) ) = 27TT S E E E ML'm'MS'm" E (a^t^wfw^') . (37) 

m— — l k .k' — W 1 va' l" .m" i,i' — l 

We note that the two set of weights corresponding to the two cleaned maps which are cross correlated are not strictly 
identical. Therefore, in the above equation, we have used a prime to distingush between the weights for two cleaned 
maps. 

We note that the primary foreground contamination at large multipole region comes from (extra-galactic) point 
sources. Diffuse galactic foregrounds are subdominant at large /. The beam deconvolution process leads to dominance 
of the effective noise of the DA maps over the unresolved point source contamination. As a result weights become 
entirely determined by the noise level of the maps and asymptotically (I — > oo) become constant for all realizations. 
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Moreover, small differeneces in weights from realization to realization due to fluctuations in noise from the mean 
level is not a concern. These fluctuations could be further suppressed by computing binned estimate of the power 
spectrum. Using simulations of the cleaning procedure with realistic model of point sources and detector noise we 
have verified that the weights remain effectively unchanged whether we use point sources or not. Since the pairs of 
cleaned maps that are cross-correlated have uncorrelated noise, the corresponding weights w kl and w' k 1 could be 
treated as uncorrelated with one another. Also, weights are effectively uncorrelated with unresolved point sources 
as they are determined by the detector noise only. In this case the residual unresolved point source contamination 
becomes 

m—l k,k ; —9 i,i'=4 

( C ^) = 27TT E E E E M L' m >M&„ m „ £ (af m ,a$«,) • (38) 

m— — lk,k' — ll',m'l",m" — l 

Next we assume that the distribution of unresolved point sources is statistically isotropic over the full-sky i.e. 

Here Cf, s ^ is the point source model as supplied by the WMAP team. This simplifies the expression for the residual 
unresolved point source contribution, 

m=l fe,fc'=9 z,i'=4 

cfSii ) - 2irr E E E Mf ml , m ,MSm' E c l s(ll ' ] ■ (40) 

m— — l k,k' — l I' ,m f i,i'=l 

The residual unresolved point source contamination then could be written as, 

m—l k,k —9 

^' S(ij) ) = 27TT E E E M lml'm' M lml'm' (w*C{fW*') . (41) 
m— — l k,k' — l V ,1m' 

Here, W} 4 is the row vector of weights for sky part k as in eq. ([7]). The elements of the point source power spectrum 

matrix are given by Cf,^' . We note in passing that this matrix is not symmetric because we are interested in 
residual point source power in the cross combination of two maps obtained from two cleaned maps which are linear 
combinations of (K, Q, V, W) and (KA,Q,V, W) respectively. Explicitly we use the following form of the matrix 



,ps(ij) _ 



/ (jKKA Q K Q (jKV qKW \ 

qKAQ qQQ (jQV (jQW 

(jKAV (jQV (jVV (jVW 

(jKAW (jWQ (jWV (jWW J 



(42) 



We note that, in general, point sources may not be perfectly correlated in all the frequencies. The point source power 
spectrum assumed by the WMAP team and used in this work assume the existence of a single point source template 
which may not be perfectly true. However, our point source correction method could easily incorporate the extra 
information of point source decoherence from frequency to frequency in the matrix in eq. (|42[) . 

From a detailed study of a single iteration cleaning method we verified that a residual point source bias in the 
auto or cross power spectrum of the cleaned maps could be very well approximated by WC[ 5S Wj r , without a need 

to compute (wC^wA A detailed discussion of this analysis and the corresponding simulations are given in 



appendix [Dl Thus we propose (Cf ) « Cf s . This simplifies the expression for the residual unresolved point 



source contribution 



=1 k,k'=9 



2/ 



l — E E E M ^'-' M ^ m 'W^c? s wjr. (43) 



m——l k,k'—l V ,m' 



Eq. [531 further be written as 

/c— 9 / m—l t 

C! Sm = E E MtuW^wf + ^— £ J2 E M? ml , m MS m W$Cf,*W>f T ] . (44) 

k=l I' \m=-lk,k'(k^k')l'.m' 
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Here Mw is the mode mode coupling matrix. The first term of the bracket is obtained in analogous method as shown 
in Ref. [34]]. The last term can also be further simplified. For this purpose we note that 



M, 



Iml' 



E 



l"m" 



(2Z + 1)(2Z' + 1)(2Z"+1) 



4tt 




I V V' 



(45) 



where w\ m are the spherical harmonic coefficients from the k th mask. Hence 

E <ml> m >MS m , = ^ + + E E ^m««&»™ ((2'" + !)(«"' + 1)) 



I" m" V' 




l V V- 



E 



i r z" 



i r i" 



m -m' m" 

mm' \ 



(46) 



Using the property of the Wigner 3jm symbol we have 

(2f + l)(2f' + l) 



EM k M* k ' 



4/T 



E E 4w-4« m «'((2i"+i)(2r"+i)) 

l"m" l"'m"' 



Hence 



2Z- 



1 ^ 

mm' 



Iml' m' Iml' m' 



(2i' + l) 



4-7T 



E E 

l"m" l"'m"' 



/ z z' r \ . . 

X dl"l'"dr, 

\ / 



1/2 / i i' z" 
^ooo 

i 



21" + 1 



'j v , m ,,w v . 




This is easily written in terms of the cross power spectra wfn of two masks k, k' 



1 V I" 




We have defined Mfi k as the cross coupling matrix of two masks k,k' . Finally we can write eq. (|3"T|) as 



cr (,j) =EE M «'Wpc{rw 



fc=l !' 



k,k'k^k' 



E^'w^c^wjK 1 



(47) 



(48) 



(49) 



(50) 



which estimates unresolved residual point source contamination in individual cross power spectrum. After correcting 
each of the cross power spectra we form a simple average to obtain final point source corrected spectrum. We note 
that the second term on the right hand side of the above equation gives entirely negligible contribution to the total 
estimate of the residual point source correction. This is a result of the fact that weights are effectively determined 
by the detector noise in the large I limit. Two different sky parts have uncorrelated noises. So weights from two 
different sky regions are uncorrelated at large I limit. In Ref. (2^ |. we have shown that the residual point source 
contamination is significantly smaller than the contamination arising from K, Ka, Q or V band. In fact point source 
residual is intermediate between V and W band point source power. Figure [5] shows, residual unresolved point source 
contamination Ylk=i Sz' Mjj,Wfi l Cf, for different values of Here i,j are the index representing the 4 DA. 

All possible combinations are explicitly shown in eq. (|42p. The total unresolved point source spectrum is shown 
as the pink line with filled circular points. The dominant contributors to the total unresolved point source spectrum 
at large / are the WW, VW and VV combinations. This is expected since the V and W band share most of the 
weights at large I because of their higher angular resolutions. Weights are negligible for K, Ka, Q bands in the large 
I limit and hence point source contamination from these bands is heavily suppressed. 

The basis of the WMAP team's 1 year power spectrum are the 28 cross power spectra which are available from the 
LAMBDA website in the unbinned form. These 28 cross power spectra are not corrected for the residual unresolved 
point sources. Following the WMAP team's 1 year bins we compare these 28 cross power spectra in Fig. [6] with 24 
cross spectra obtained from our own 1 year results. These 24 cross spectra also are not corrected for the residual 
unresolved point sources and show very little dispersion compared to the WMAP results. The 'uniform average' power 
spectrum plotted in green line in fig. [7] has less excess power near the second acoustic peak compared to an 'uniform 
average' of the WMAP's 28 cross spectra (red line). This merely shows that we have removed some amount of point 
sources at this range of / during our cleaning. 
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FIG. 7: Comparison of the WMAP's average power spectrum (red) with our binned power spectrum without any point source 
subtraction ( green). Clearly starting from the first acoustic peak we have less point source contamination. WMAP's final 
binned power spectrum is also shown in blue for comparison. Interestingly the notch at 1=4 appears to be reduced in the 
WMAP's average power spectrum (red) . 
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D. Computing error bars 

We rely upon Monte-Carlo simulation to compute the error bars on our power spectrum. We generate synchrotron, 
free-free and thermal dust maps corresponding to different frequencies in a given combination using the publicly 
available Planck Sky Model. Each of the random realizations of CMB and foreground maps are convolved by the 
appropriate beam function for each detector. Random noise maps corresponding to each detector are generated by 
first sampling a Gaussian distribution with unit variance. In the final step we multiply each Gaussian variable by the 
number ctq/ a/ N p to form realistic detector noise maps. Here <7o is the noise per observation of the detector under 
consideration and N p the effective number of observations at each pixel. These realistic maps with detector noise, 
foreground and CMB signal are then passed through the cleaning pipeline. The error-bars for our power spectrum 
correspond to the standard deviation of the power spectrum obtained from Monte-Carlo simulations . 

Due to the (Kp2) mask applied to remove potential foreground contaminated regions, the neighboring multipoles 
become coupled. In the presence of detector noise, correlation between neighboring Cf becomes stronger. The 

covariance matrix ^ACfACf, ^ = ((Cf — (<Of^)(C£, — (Cf,\)^ obtained from simulations is therefore expected to 
have non-diagonal elements. It is convenient to bin the power spectrum in order to minimize the correlations and 
errors. We have considered a binning identical to that used by the WMAP team in their analysis. Let Cb denotes 
the binned power spectrum. Then the covariance matrix of the binned spectrum is obtained as (ACfACf, 



(C§ — {C^j){C^t — (Cyj)j. The standard deviation obtained from the diagonal elements of the binned covariance 

matrix was used as the error-bars on the binned final spectrum extracted from the WMAP data. We also define a 
normalized covariance matrix of the binned power spectrum following, 

(AC b ACv) 

C bb ' = j N 1 , (51) 



{(AC b y) ((AC* 



I - 



wherein all the elements of this matrix are bound to lie between [—1, 1]. This correlation matrix represents the actual 
bin to bin correlation matrix following the cleaning method. In the left panel of the figure [5] we show the correlation 
matrix for WMAP 1 year simulations. The right hand panel of this figure is the corresponding plot for the WMAP 3 
year analysis. Both these matrices are seen to be sufficiently diagonal dominated. 



IV. RESULTS 



Figure [H] shows the main result of the power spectrum for CMB anisotropy estimated using our analysis of WMAP 
1 year and WMAP 3 year data. The blue line shows the WMAP 3 year power spectrum and the red line shows the 1 
year spectrum. All these spectra are corrected for residual unresolved point source contamination. In the lower panel 
of this figure we show the residual unresolved point source contamination for 1 year and 3 year respectively. In what 
follows we describe the power spectrum obtained by us for WMAP 1 year data and WMAP 3 year data respectively. 



A. WMAP 1 year data 



Using the 48 foreground cleaned maps obtained from the WMAP 1 year maps we obtain a 'Uniform average' power 
spectrum following the method mentioned in section IIII Bl The residual unresolved point source contamination is 
removed following IIII CI The estimated power spectrum with error bars is plotted in red in figure [9] 

We find a suppression of power in the quadrupole and octupole moments consistent with the results published by 
the WMAP team. However, our quadrupole moment (146/iif 2 ) is little larger than the quadrupole moment estimated 
by WMAP team (123^if 2 ) and Octupole (Abb^iK 2 ) is less than the WMAP team result (611^if 2 ). The 'Uniform 
average' power spectrum does not show the 'bite' like feature present at the first acoustic peak in the power spectrum 
reported by WMAP [Io| . We perform a quadratic fit of the form AT; = AT; + a(l — l n ) 2 to the peaks and troughs 
of the binned spectrum similar to WMAP analysis [39j . For the residual point source corrected power spectrum we 
obtain the first acoustic peak at I — 219.8 ± 0.8 with the peak amplitude ATi = 74.1 ± 0.3{iK, the second acoustic 
peak at / = 544± 17 with the peak amplitude AT; = 48.3 ± \.2\iK and the first trough at / = 419.2 ± 5.6fiK with peak 
amplitude AT; = 41.7 ± lfiK. The left panel of figure ITD1 shows the three different ranges of multiples used to find 
out peak and trough positions and their corresponding amplitudes AT;. (A similar plot for the WMAP 3 year data 
is shown in the right panel of this figure. The results for WMAP 3 year analysis are summarized in section TlV Bl ) 




As a cross check of the method, we have carried out the analysis with other possible combinations of the DA maps. 

1. The WMAP team also provide foreground cleaned maps corresponding to Ql to W4 DA (LAMBDA). The 
Galactic foreground signal, consisting of synchrotron, free-free, and dust emission, was removed using the 3- 
band, 5-parameter template fitting method [ill ]. We also include K and Ka band maps which are not foreground 
cleaned. The resulting power spectrum from our analysis matches closely to the 'Uniform average' power 
spectrum. 

2. Excluding the K and Ka band from our analysis we get a power spectrum close to the 'Uniform average' results. 
Notably, we find a more prominent notch at I = 4 similar to WMAP's published results. 

In case of 'Uniform average' a maximum difference of 92 fj,K 2 is observed only for octupole. For the large multipole 
range the difference is small and for I — 752 it is approximately 48 /j,K 2 . This is well within the la error bar (510/iX 2 ) 
obtained from the simulation. This shows that our foreground cleaning is comparable in efficiency to that obtained 
by employing template fitting methods that rely on a model of foreground emission to estimate the contamination at 
the CMB dominated frequencies. 




FIG. 9: Comparison of 1 year 4 channel with that of 3 year 4 channel power spectrum. The best fit WMAP's power spectrum 
is shown in black line along with cosmic variance band. The bottom panel of this figure shows the residual unresolved point 
source contamination for both the power spectra. 



The Monte Carlo simulations of our cleaning method also reveals the negative bias in the low I moments. The origin 
of this negative bias is explained in section Hi Bl For I = 2 and I = 3 the bias is —27.4% and —13.8% respectively. 
However this bias become negligible at higher I, e.g. at I = 22, it is only —0.8%. This bias can be explained in terms 
of an anti-correlation of the CMB with the residual foregrounds in the cleaned map. For further details of the bias 
we refer to appendix [El The standard deviation obtained from the diagonal elements of the covariance matrix is used 
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FIG. 10: The left panel shows the 3 different multipole ranges used to obtain positions of the first peak, the first trough and 
the second peak from our point source subtracted power spectrum using 1 year WMAP data. Before fitting the 1 year power 
spectrum was binned in the same manner as the WMAP's binning. The box error-bars are used to indicate x and y error-bars. 
The right panel shows the same figure but using the WMAP 3 year data. A correction due to residual unresolved point sources 
was performed prior to fitting. 





FIG. 11: The left panel shows (in red points) ensemble averaged power spectrum from 110 Monte Carlo simulations of our 
power spectrum estimation method. The simulations were carried out using the 1 year WMAP detector noise maps available 
from the LAMBDA website. We use publicly available Planck Sky Model to generate the diffuse foreground models. The 
recovered spectrum is binned in the same manner as WMAP 1 year power spectrum. The input theoretical spectrum is shown 
in blue line with cosmic variance. The right panel is same figure but with 3 year noise maps. The 3 year noise maps were 
generated following the method described in the text. The spectrum is binned following the binning scheme of the WMAP 
team's analysis of 3 year data. 



as the error bars on the CVs obtained from the data. The ensemble average of 110 cleaned power spectrum is shown 
in the left panel of the figure [TT1 The presence of bias at the low multipole moments is visible is this figure. 



B. WMAP 3 year data 

1. 4 channel combinations 

We analyse the 3 year WMAP data by using a procedure identical to that used for the 1 year data. The details 
are given in section [TTJ In figure [121 we show each of the 24 cross power. The 'Uniform average' of these 24 cross 
power spectra is also shown in this figure in red line with blue error-bar. This figure is similar to the 24 individual 
cross power spectra obtained for WMAP 1 year data [U, [22[ ■ For the 3 year data all the 24 cross power spectra show 
very little dispersion till the second trough. In comparison the 1 year data shows small dispersion only till the second 
acoustic peak [2lJ. This may be explained due to the effectively lower detector noise in the 3 year data compared to 



20 



the 1 year data. 
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FIG. 12: The 24 cross power spectra for the 3 year WMAP data are shown with the detector noise bias removed. The red 
line with blue error-bars is the 'Uniform average' power spectrum. In the inset of this figure we show all the possible cross 
correlations of the cleaned maps that give rise to the 24 cross spectra. 

We perform a peak fitting to the peaks and troughs of the 3 year power spectrum as well. A parabolic function 
of the form AT; = ATj + a(l — Iq) 2 is fitted to the peaks and troughs of the power spectrum amplitude. We use 
the 3 year binned data for the fitting purpose. For the point source corrected power spectra, the first acoustic peak 
has an amplitude AT; = 74.4 ± 0.3/i-K" at the multipole position I = 219.9 ± 0.8. The first trough is located at 
I = 417.7 ± 3.2 with an amplitude ATJ = 41.4 ± 0.6/iK. The amplitude and position of the second acoustic peak are 
given by AT = 49.4 ± 0.4/iif, I = 539.5 ± 3.7. In the right panel of figure fTUl we show the different multipole ranges 
used to obtain the position of the peaks and troughs and their amplitudes. 

2. 3 channel combinations 

We also follow an alternative approach in which we use only the Q,V and W band DA maps. This is similar to 
the 4 channel cleaning, however now we only get 24 cleaned maps. We form a total of 12 cross power spectra from 
these cleaned maps after applying 3 year Kp2 mask. After debiasing all the power spectra by the coupling matrix 
and removing beam and pixel effect we obtain an '3 channel uniform average' power spectrum. We find that the 3 
channel spectrum matches well with the 4 channel spectrum. 

C. Comparison of 1 year and 3 year power spectra 

Fig. [9] compares the power spectra obtained from the 1 year and 3 WMAP data binned identically using the 
WMAP team's 1 year binning method. They match closely with one another and also with the WMAP's best fit 
power spectrum available from the LAMBDA website. For the residual point source corrected power spectrum we 
obtain for WMAP 1 year (WMAP 3 year) data the first acoustic peak at I = 219.8 ± 0.8 (219.9 ± 0.8) with the peak 
amplitude AT; = 74.1 ± 0.3^ (74.4 ± 0.3fj,K), the second acoustic peak at I = 544 ± 17 (539.5 ± 3.7) with the peak 
amplitude AT; = 48.3 ± \.2\iK (49.4 ± OA^iK) and the first trough at / = 419.2 ± 5.6/iif (417.7 ± 3.2) with peak 
amplitude AT; = 41.7 ± 1 fiK(UA ± 0.6fiK). 

We note that our cleaning method significantly removes unresolved point source contamination. The original I 
dependence of the unresolved point source power spectrum present in the foreground contaminated maps (as well 
present in template cleaned maps) is significantly reduced and becomes independent of I at large I. We also note that 
the unresolved residual point source contamination is less by about w 50[iK 2 in 3 year power spectrum than the 1 
year power spectrum. This is expected. The WMAP supplied 3 year Kp2 mask removes more point sources than the 
1 year Kp2 mask. 
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V. CONCLUSION 



The rapid improvement in the sensitivity and resolution of the CMB experiments has posed increasingly stringent 
requirements on the level of separation and removal of the foreground contaminants. We carry out an estimation of 
the CMB power spectrum from the WMAP data that is independent of foreground model. The method does not rely 
upon any foreground template and employs the lack of noise correlation between independent channels. This paper 
is a detailed description of the first estimate of the CMB angular power spectrum solely based upon the WMAP 
data. In this paper, we present an indepth study of the biases that arise in the foreground cleaning. In particular, we 
provide an understanding and correction for the negative bias at low multipoles reported in our earlier work [2lj . 

Usual approaches to foreground removal, usually incorporate the extra information about the foregrounds available 
at other frequencies, the spatial structure and distribution in constructing a foreground template at the frequencies of 
the CMB measurements. These approaches may be susceptible to uncertainties and inadequacies of modeling involved 
in extrapolating from the frequency of observation to CMB observations that a blind approach, such as presented 
here, evades. The understanding of polarized foreground for CMB polarization maps is rather scarce. Hence modeling 
uncertainties could dominate the systematics error budget of conventional foreground cleaning. The blind approach 
extended to estimating polarization spectra after cleaning CMB polarization maps could prove to be particularly 
advantageous. 
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APPENDIX A: ANALYTIC DERIVATION OF WEIGHTS AND CLEANED POWER SPECTRUM 

The main idea behind the blind foreground cleaning method used here is entirely based upon the minimization of 
total power in the cleaned map in multipole space [18|. Weights for different channels are obtained minimizing the 
total power 

(jciean f t h e cleaned map. However, we ensure that the CMB angular power spectrum is conserved during 
cleaning by imposing the constraint Wieo = W T = 1 on the weights. The solution for the weights that satisfy these 
conditions is the point in weight space where normals to the functions /(Wj) = WiCiW^ and .g(Wi) = Wieo are 
parallel to one another. Following Lagrange's multiplier method, this is cast to an equivalent problem of minimizing 

WiCiW^ - AWie . (Al) 

Here A is the unknown Lagrange multiplier parameter which could be determined from variational principle. At the 
extrema, the expression in eq. (|A1[) is unchanged under small variations in Wj leading to 

AWiCiW^ + WiCi AW^ - AAWie = . (A2) 

Since the power spectrum matrix, Ci is a symmetric matrix, the first two terms of the left hand side are equal to one 
another. Hence we obtain, 



AW 



2CiW? - Ae 



. (A3) 



8 The HEALPix distribution is publicly available from the website http://healpix.jpl.nasa.gov 
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Since this relation is true for any arbitrary variation AWi, we obtain 

2CiW! T - Ae = . (A4) 

We introduce a (non zero) square matrix Gi. Later we will identify Gi as the Moore Penrose Generalized Inverse 
(MPGI) of the covariance matrix Ci. After multiplication from left by this matrix we can rewrite the above equation 
as 

2GiCiW 1 T - AGie = . (A5) 

Hence we obtain 

A= 2ejG^WT (A6) 
e Gie 



Now we use the constraint Wieo = 1. Assuming Wi ^Owe multiply eq. (IA4p from left by Wj and to obtain 

2WiCiW^ = AWie = A . (A7) 

Using eq. (|A6[) and eq. (|A7|) we obtain 

2W lCl WF = 2 e °^'^' , (A8) 
e^Gie 

or, 

(w.-^g^QW^O. (A9) 

If we neglect solutions which belongs to null space of Ci, 9 i.e., assuming CiW^ ^ 0, we obtain 

W. = ■ (MO) 

e^Crie 

We impose the restriction on Gi that it is symmetric (since the MPGI of a symmetric matrix is symmetric) to obtain 

W? = 4^- (AH) 
The corresponding power spectrum of the cleaned map is given by 

W lCl W^ = _^L Cl ^- = — i— , (A12) 
ejGie ejGie ejGie 

where we have imposed the condition GiCi = CjGi. It is easy to note that the choice of symmetric Gi helps to 
obtain a simplified expression for the cleaned power spectrum. Also one can verify that, if G1C1G1 = Gi is satisfied 
and Gi is symmetric then Gi satisfies all the defining conditions of MPGI of Ci. Hence, we can identify Gi as the 
Moore-Penrose Generalized Generalized Inverse (MPGI) of Ci (and vice versa). 

Eq. (|A12|) has an interesting property. It remains valid even when the full covariance matrix Gi is singular. A 
singular full covariance matrix is encountered in noiseless case (or numerically very low noise case) when the foreground 
components follow a rigid frequency scaling and total number of components (all foregrounds and CMB) become less 
than the number of available channels. If Ci is non singular we can replace Gi everywhere by Cj -1 . This is because 
MPGI of a non-singular matrix is its inverse. 



There is a physical justification behind neglecting this solution. The weights satisfying CiWj 1 " = do not preserve CMB power while 
minimizing total power in the cleaned map. These solution for the weights merely sets the total power in the cleaned map WjCiW^ = 0. 
Thus we are not interested in the space of solutions to CiW^ = 0. 
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APPENDIX B: PARTITIONING THE SKY 



An important advantage of the foreground cleaning analyzed here is that we can allow the weights to vary with 
sky positions as well as with the multipole moment. To allow the weights to vary with sky positions we can partition 
the sky into several regions depending upon the level of foreground contamination. (Alternatively the partition could 
be done directly using the knowledge of expected spectral index dependence on the sky.) We followed the procedure 
of Ref. [IH to partition the sky. Each of these partitions are identified with a sky masks. The mask takes non zero 
value at all pixels contained within the sky partition represented and is zero outside. In this section we describe the 
procedure of constructing these masks. 

There are 2 Difference Assemblies for Q band, 2 for V band and 4 for W band. To make the masks we first averaged 
all the DA maps for a given frequency band. Correspondingly we averaged the beam functions for each frequency band. 
For the K and KA bands there are only one difference assembly in each case. Therefore for these bands no averaging 
was done. We smoothed all the five maps (e.g. K, Ka, Q, V, W) by the resolution function of the K band, which 
has the lowest resolution. We obtained four difference maps W-V, V-Q, Q-K, K-Ka, out of these 5 smoothed maps. 
The smoothing was performed first obtaining the a; m coefficients. As each averaged map was effectively smoothed by 
the averaged beam function (corresponding to each channel) during the observations, we decided to first remove the 
beam effect. Then we smoothed each map by the common resolution of the lowest frequency band. Mathematically, 
the difference maps were obtained as follows: 

-pK tjK 

„W-V W a l V a l /ni *\ 

a hn - a lm-^W ~ a lm-^v~ ' V C 1 ) 



V-Q 
Hra 



a 



v 



Im r/V 



Hra 



B? 



(B2) 



„Q-K — „Q £iL_ _ „k 

"I 



(B3) 



K — KA _ K „KA °l / nA \ 
a lm ~ a lm ~ a lm ~^KA ' i B4 -> 

B l 

The ai m were converted to difference maps using HEALPix supplied subroutine alm2map. 

Next we construct a junk map out of these four difference maps assigning at each pixel the absolute maximum 
value among the four difference maps. We down-sample the junk map using the HEALPix supplied program udgrade 
to a resolution of N s id. e — 64. We identify 7 different sky mask partitions from this low resolution junk map after 
applying cutoff corresponding to the following temperature thresholds (in fiK) T > 30000, 30000 > T > 10000, 
10000 > T > 3000, 3000 > T > 1000, 1000 > T > 300, 300 > T > 100 and T < 100. The partition T > 30000 
is maximally contaminated by foreground emission. The second dirtiest partition is disjoint on the sky and we use 
3 separate masks corresponding this partition. The resulting 9 masks are then converted back to the HEALPix 
resolution N s n e = 512 using udgrade routine of HEALPix. Next we smooth each mask using a Gaussian beam of 
FWHM 30' and redefine smoother mask boundaries at the threshold of 0.5. 10 In this case we found that almost the 
entire sky is covered by the 9 masks (except for a few pixels in the sky). 

To convert the second dirtiest region into 3 mask files we upgraded the -/V s id e = 64 resolution map to the -/V s id e = 512. 
Then using IDL task mollcursor, we found out the extension in galactic 8, <f> coordinates of the 3 different parts of the 



We found that if we smooth them by a Gaussian function of FWHM 2 degree (as mentioned in Ref. |18|) and then define the boundaries 
at the threshold of 0.5 the resulting 9 masks cannot cover the entire sky. In that case there are some regions near the galactic plane 
which do not belong to any of our 9 masks. Obviously if such regions are not covered by any of the masks and remain present in our 
final cleaned map then final power spectrum will be contaminated by the foreground. Therefore we chose to use a Gaussian function 
of lower FWHM of 30' for smoothing. The reason why in our case, masks smoothed by a 2 degree Gaussian function cannot cover the 
entire sky is clear. This is actually dependent on the common beam function by which we are smoothing each map before forming the 
difference maps. We found that the masks near the galactic regions contain a few small isolated regions. Therefore smoothing by a 
Gaussian with FWHM as large as 2 degree, gives this isolated small regions maximum value far less than unity, in fact maximum value 
becomes quite near 0.5, the cutoff value, after smoothing. Consequently applying a cut of 0.5 removes most of the part of these isolated 
regions, leading to some part of the sky, uncovered by the masks. 
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second dirtiest region. We then implemented a method which determines the pixel index in 'ring' format for these 3 
regions. Finally we converted them to 3 different masks at N B ^ B = 512. The masks were then smoothed by Gaussian 
function of FWHM =30'. We applied a cutoff of 0.5 to each of them. The 8 th , 7 th and 6 th masks were numbered 
according to the descending order of maximum pixel value in the junk map at resolution A^dc = 512. 

Figure [2] shows our 9 different masks that partition the sky based on estimated level of foreground contamination. 
These regions are similar to what is shown in Ref. [HI]. As the figure shows, one side of the band near the galactic 
plane is more severely foreground contaminated. 



APPENDIX C: COMBINING CROSS POWER SPECTRA 



The basis of our final power spectrum are a set of 24 cross power spectra where the detector noise bias has been 
removed. An uniform weighting of cross spectra is used to obtain the final power spectrum. In this section we describe 
the procedure to obtain the 24 cross spectra and their combination to form the final spectrum. 

Let us assume that AT l (9, (f>) represents one of the 48 final cleaned map, where (i = 1, 2, 3, 48). To remove the 
residual foreground contamination near the galactic plane (figure [3]) we apply the Kp2 mask supplied by the WMAP 
team on each of these 48 maps. The masked map can be represented as 

AT'\n) = W(h)AT\h). (CI) 

Here W{6, 4>) represents the Kp2 mask. In the next step, we obtain a cross power spectrum between pairs of foreground 
cleaned maps that have uncorrelated noise (recall that noise in different DAs are uncorrelated) . If a\ m and a\ m are 
the spherical harmonic coefficients obtained from two such maps 



W(n)AT ij '(n)dn. (C2) 



The cross power spectrum is obtained using, 



■m—l 



7, J 4 



C]' = > " . (C3) 

1 ^ 21 + 1 y ' 

m= — I 

Here the super scrip t ij represents the cross power spectrum obtained from the i th and j th foreground cleaned maps. 
Following Ref. [3J], the ensemble average of the C/ estimated from the partial sky is related to the C; from the full 
sky as 

C;) = ^M„,((7f,} , (C4) 
v 

where Mw is a coupling matrix. This matrix represents the fact that when we multiply our map by the weight 
function (Kp2 mask) in the pixel space we are effectively performing a smoothing operation and neighboring spherical 
harmonic coefficients get coupled. An analytic expression for this coupling matrix is given in Ref. [341 ] . 

ls—\h— h\ \ / 

where the last term is the Wigner-3j symbol and Wi is the power spectrum of the mask under consideration. Although 
in our case we are interested in the partial sky power spectrum which has been obtained by cross correlating two 
different maps, it is easy to show that in this case also eq. (|C4|) remains valid. The rank of the coupling matrix is 
limited to l m ax = 1024. We use numerical routines from the NETLIB package [4(| to compute the Wigner -3j symbol. 

Eq. (|C4p is true only for an ideal observation with infinite angular resolution. In practice all instruments have finite 
angular resolution given by the beam function B(9 1 </>) of the instrument. Mathematically, we may write the observed 
temperature anisotropy, 

AT(n) = J AT'{n')B(n,h')dn n , , (C6) 

where AT'(n') is now the full sky map in the absence of the beam. In most of the CMB experiments the beam 
function is circularly symmetric to a good approximation, i.e., depends only on the angle 6" = cos _1 (n • h') between 

B(n, h') = B(n ■ ti) . (C7) 
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We may expand this function in terms of Legendre polynomials, 

B(n-ti) = J2^T— B i p i(n-n')- (C8) 



47T 

1=0 



Here, Bi's are Legendre transform of the beam function. Substituting this in eq. (|C6[) and expanding AT and AT' in 
spherical harmonics, we obtain 



Z— oo m— I I — oo m — „ I' 



I, IJfJ lib 1 L UXJ III 1 p I IAJ C17// I 

£ aim^m(n)= 51 £ a 'l'm> / ^'m'(0 E ~ ^~ B '" P '" (" ' • (C9) 

Z=0 m=-i i'=0 m'=-C Z"=0 

Using the addition formula 

4tt m "= l " 

P v^-ti) = Wn^l ^2 Y v/rnt/ (h)Yfi mt ,(h') (CIO) 

m"= — I" 

and the orthonormality property of spherical harmonics 

f yjw(W»m»(«') ='V/-\,,',,,- (01.1) 

we obtain, 

l—oorn—l l' — oo m' — V 

J2J2 a ^ Y ^=J2 E CWPi-. (Ci2) 

Z— m=— Z m r =—V 

We again use orthonormality of spherical harmonics to obtain 

a; m = aJ m .-B ; • (013) 

This relation shows that due to the finite resolution of the instrument spherical harmonic coefficients get multiplied 
by the Legendre transform of the beam function. Hence it is easier to account for the effect of a circular beam in the 
spherical harmonic space than deconvolving the map by the beam function in the pixel space. 

The effect of finite pixel size of the map has similar effect on the recovered spherical harmonic coefficients. The 
pixel window functions pi have been supplied with the HEALPix distribution for different resolutions. Taking into 
account finite pixel size we have 

a lrn = a'l m BlPl . (C14) 

Hence the recovered power spectrum is related to the actual CMB sky power spectrum as 

Ci = ClBfpf . (C15) 
In presence of both beam and finite pixel effects, eq. (|C4[) obtained for the ideal case 11 is now modified to 

[Ci) =Y, M w {Ci) B*pI . (C16) 
I' 

A generalization of this expression in case of cross power spectrum has been also reported in Ref. [35| . The partial 
sky cross power spectrum is related to the full sky power spectrum in the following manner 

C\t)=Y,M lv (ct)B\,Bi,pl. (C17) 



11 In Ref. pS4l | the authors followed the convention that Bi represents combined smoothing effect due to pixel as well as beam. Our notation 
is different from them. 
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Here we note that the noise terms drop out in the cross-correlation procedure. 

The final estimate of the full-sky spectrum is obtained by simply inverting the coupling matrix. We have used a 
singular value decomposition technique to invert our matrix. We checked that this matrix gives unbiased estimates 
of the full sky power spectrum using Monte Carlo simulations of CMB maps and using Kp2 cut. 

The final full sky estimate of the power spectrum is obtained from 

C\ = Y,{M- l ) lv Cii/BtB{pl (C18) 
v 

with i = i = (1,2,3,4, ...,24) as there are 24 cross-correlations possible that do not have detector noise bias. 

Using the 24 cross-power spectra we form an averaged power spectrum following 

i=24 

Q = J2 CjN cross (C19) 
where N cross = 1/24. Now we bin them in the same manner as the WMAP. The binned power spectrum is defined as 

with Al = l m ax + 1 — l-min- Eqn. (|C20j) defines our final power spectrum. 



APPENDIX D: ESTIMATION OF RESIDUAL UNRESOLVED POINT SOURCE CONTAMINATION 

In this appendix we study the point source bias in the auto power spectrum of a cleaned map and in the cross 
power of two cleaned maps. From a detailed study of a single iteration cleaning method we verify that a residual point 
source bias in the auto or cross power spectrum of the cleaned maps could be very well approximated by WC[ )s Wj r , 

without a need to compute ^WCpW^. For the auto power spectrum we optain an analytic point source bias in 

terms of the noise covariance matrix and the point source covariance matrix. We note that, in practice when we 
estimate cosmological power spectrum from the WMAP data we do not use the auto power spectrum of a cleaned 
map. Nevertheless the auto power spectrum is important to study along with the cross power spectrum for a deeper 
understanding of the point source bias. We note that in all the Monte Carlo simulations in this section we treat point 
sources as fixed templates. The sources of randomness come from CMB and detector noises. 

1. Auto Power Spectrum 

Following eq. (|27p the cleaned auto power spectrum obeys, 



c Clean ) = (df) + ( j— ) + (1 - Tl c ) . ( 1) 1 ) 

Clearly the positive bias is given by ^l/eJ(C[ +N ) eo^} which is caused by the detector noise and foreground 

covariance matrix of the map. In this appendix we are basically interested in studying point source bias and since 
diffuse foreground is not a concern at large I, we do not include any diffuse foreground component in this study. 
Considering the second term on the right hand side of eq. (|D1|) we note that it is not possible to identify a bias 
which solely comes from point sources. Another point to note is that we know only the mean noise covariance matrix, 
, not the empirical noise covariance matrix, Cj , where Cf = Cf + 5Cf and SCf 1 denotes noise fluctuation on 

the true noise level. Therefore we express (VeJ(Cf +N ) 1 e ) in terms of a mean noise covariance matrix and the 

foreground covariance matrix. As we will see such a simplification helps us to obtain an analytical expression for the 
point source bias. Keeping in mind that WMAP detector noise level is much larger than the point source spectrum, 
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we simplify yl/ej (C[ +N ) eoy analytically by making expanding to first order 12 , in terms of C[ >S (C{ SI ) , 

(C[+ N ) _1 = (C>f * - (CH'^rCCf)" 1 , (D2) 

Due to the low point source contamination relative to the WMAP noise level it is reasonable to assume that 
(eJ(CN) _1 CP s (Cf I ) _1 eo)/(eJ(Cf , ) _1 e ) < 1. Expanding to first order we obtain, 

/ l\ + / - T «yr' cr ) , (D3) 

\eJ(C[+ N )- 1 e„/ W(C~f W W(Cff e„ eJ(Cf ) e„ / 

We interprete the first term on the right hand side as the detector noise bias. The second term is the point source 
bias. In the following subsections we recast these results in terms of the mean noise covariance matrix. We also verify 
our analytic expressions using Monte Carlo simulations. 

a. Noise induced bias 

We first turn our attention to the detector noise bias. We estimate this term upto second order in (JC^C™ 



(e?(CN)- 1 ,5CN(CN)- 1 eo) ; 
(eJ(CN)- 1 e ) 2 



(D4) 



The noise fluctuations are assumed to have zero mean, /^Cj^y = 0. Hence the second term in the bracket vanishes 

on the ensemble average. Only the numerators of the third and fourth terms in the bracket on the right hand side of 
eq. (|D4[) are stochastic variables. On the ensemble average the numerator of the third term becomes, 



/ 

\i,j,k I 

where we have assumed that the noise covariance matrix of WMAP is diagonal. We simplify the ensemble averaged 

quantities using the relation eq. (28) of Ref. 35]. If we consider a general term of the form ^SC^^ 5C^ kp ^\ then 

for uncorrelated detector noises, the ensemble averaged quantities will survive only for the same pair, i.e. when 
= (k,p). We further note that, 

{ 21+1^1 °Z ' 1 — J 

Using above relations the third and fourth term in the bracket of eq. (|D4[) become 

/eT( C N)- 1 <5Cr'(CN)- 1 ( 5Cf I (Cf I )- 1 eo\ _ (l + n c ) 



3 J(Cf I )- 1 eo / 2/ + 1 



(D6) 



Following a similar method we also simplify the fourth term in eq. (|D4|) . The final result is 

/ (e?(CP)- 1 ^(Cr J )- 1 e ) 2 \ 2 
\ (e^CN^eo) 2 / 21+ 1' 



(D7) 



We note that we cannot use Sherman-Morrison formula to decouple the foreground and detector noise bias from the term 
— - — — 1 — ; in an useful form. 
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FIG. 13: The difference between the empirical noise bias and the theoretical estimate of noise bias given in eq. ID8l for V and 
W bands. 



Using eqs. ijDlj) . (|D6|) . (|D7|) wc obtain 




2l + 2 -^ (D8) 



\eJ(Cf)- 1 eo/ \eJ(CN)- 1 eoV 2/ + 1 2l + l)/ e^CfT'eo 2? + 1 

This enables us to compute the noise bias in terms of the theoretical noise models. We verified this equation by 
Monte-Carlo simulations of 1000 noise covariance matrices corresponding to WMAP V and W bands. The noise maps 
from which these empirical noise covariance matrices were computed were formed following a method similar to that 

mentioned in section UlI Dl In Fig. [13] wc show the difference between the empirical noise bias /l/eo(Cf) eo^ 

and the second order analytical expression for noise bias, given by the right hand side of eq. (|D8|) . The difference is 
consistent with zero. The large fluctuations at large I is caused by the beam deconvolution effect in the effective noise 
covariance matrix. 



b. Point source bias 



Here we study the point source bias. We perform Monte-Carlo simulations in a single iteration cleaning with 
CMB, point sources and detector noise having WMAP's noise level. We compute ensemble averaged cleaned power 
spectrum ( (jciean \ f rom 1000 such simulations. The validity of eq. (|D8[) allows us to compute noise bias in terms of 



the theoretical noise covariance matrix and subtract this analytical result from {C^ lean j. After correcting for noise 

bias wc find that (^cciean^ _ ^i/e^Cf 1 ) 1 eo^j ((21 + 2 — n c )/(2l + 1)) still has a residual positive bias at multipole 

range I > 400. This excess is a clear demonstration of the presence of point source bias. The brown line in Fig. [TJ] 
is the power spectrum, corrected for the noise and the negative CMB bias of the form —0^/(21 + 1). The excess is 
caused by the point source bias. We analytically compute the point source bias in terms of theoretical noise covariance 

matrix similar to eq (|D8|) . Following a first order expansion in Cj 3S (Cj sr ) and a second order expansion in terms of 



SC^Cf 1 ) 1 we obtain, 



e^CrT 1 r, ps (Cn^eo \ _ eT(Cf r 1 „p S (CfT'eo 21 , tr(Cf (Cf 1 )" 1 ) 



e^Cff'eo ' e^CNpeo/ ^(C^e * eJ(Cf)^„a + l ( 2 f + l)(ef (Cf J^ep) 
Again we verify the validity of the above expression following Monte Carlo simulations. We compute the point source 

)■ 



_ C p S vj_i / = ~uv~i / c ps vri > ^ + / (D9) 



bias from 1000 cleaning simulations following (df lean ) _Cf -Cf/(2l + 1) - (l/e£(Cf) 1 e ) ((21 + 2- n c )/(2l + 1)) 
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6000 




Multipole, I 

FIG. 14: The auto power spectrum without correcting for the residual unresolved point source bias is shown in brown line. 
The blue line shows the theoretical CMB power spectrum. 



where Cf denotes the theoretical CMB power spectrum. In Fig. [15] the green line shows the binned estimate of the 
point source bias from simulations. The error-bars plotted in this figure are valid for the mean of 1000 simulations 
of random samples. The black line with filled points shows the analytical estimate of the point source bias computed 
from the right hand side of the eq. (|D9|) . These two results match closely except at the last two points. However, 
one can easily conclude from the error-bars plotted in this figure that differences between the two lines at these two 
points are not significant. We interprete these deviations as noise fluctuations. 

Given the large noise level of WMAP maps, one might expect that ^WiCj ,s Wj r ^ would be a good approximation 

of the point source bias that is given by the second term of eq. (|D3j) . In the limit Ci — » , applicable when noise 
dominates at large I and point source bias is also significant, these two terms are identical. Using the analytical form 
of weights as in eq. (fl"4)) , we find that up to the first order in Cf s Cj _1 , 

' e^Cf)' 1 cps e T ( C N)' 1 \ Cf / e?(CN)- 1 Cr(Cr)' 1 e \ 

TfrNr 1 ^ 1 T/pNr 1 „ / 2Z + i\ pT/nNr 1 



WiC^Wj 1 ) = u p '_ x c 



,eJ(C^) eo eJ(C^) e / ^ + eJ(Cf*) e / 

+^Mcr((cfT 1 ))- (dio) 

Here, Cf is the theoretical CMB power spectrum. However, the last two terms in the above equations are negligibly 
small compared to the first term on the right hand side. From simulations we find that these two terms contribute 
less than 0.1/j.K 2 to the point source residual when multiplied by the prefactor 1(1 + l)/(27r). The reason why these 

terms become negligible is that, apart from a multiplicative first order term, C[ 5S (C{ SI ) , the first term on the right 
hand side go as Cf, whereas the last two terms goes as Cf/(2l + 1). Indeed, the noise becomes much larger than 
CMB power spectrum at large / due to beam deconvolution, making the two last two terms insignificant compared to 
the first term. In figure [TBI we show magnitude of the individual contributions arising from the second and third term 
following 1000 Monte Carlo simulations of the cleaning method. We use V and W bands in the simulations using 
realistic WMAP noise levels and residual unresolved point source model. The red line shows the difference between 
the two terms which is less than O.Ol/iK 2 or even smaller at I > 400. Hence it is justified to assume, 

/WiCrwrt - / ^(C") c ps goW) \ . (Dll) 



eJ(Cf) e eJ(Cf) e , 



In Fig. [14] we show the point source bias computed following^ WiCj W/ ) in blue line. This matches well with the 
black line with filled circle, justifying eq. (|D11[) . 

The unresolved point source bias is then (WiCP s Wj r ). Now we propose that, the point source bias could also 



J \ " l 

be taken into account following WiCf^Wj 17 ~ ^WiCj >s W J [, V The reason is that, at large I where detector noise 
dominates. The weights are entirely determined by the mean noise level. Therefore weights become approximately 
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FIG. 15: The point source bias, (WiCf^Wj ) (blue line), and the approximation WiC^Wj (red line) for the case of 
auto power spectrum of the cleaned maps. The empirical point source bias computed by using (cf let 



- Ct - C?/(2l + 

l/eg(C?)~~ e ) ((2Z + 2-n c )/(2Z + l)) is shown in green line. The error bars are computed only for the term 



1/eo (Cf 1 ) eoj without including the cosmic variance. All the data points are binned following the WMAP's 3 year binning 
method. The analytical point source bias computed following eq. (|D9|) is also shown in black line with black dots. 
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FIG. 16: The contributions arising from the third (blue line) and fourth (pinkline) terms of eq. (|D10[) and their difference. 



constant from realization to realization. Small fluctuations in weights due to noise fluctuations are not important 
compared to the magnitude of total point source bias ^WiCf^Wj 1 '). The advantage of this method is that, we can 
estimate point source bias in terms of the weights which are the results of a cleaning. In the figure [151 we show how 
well WiCj^Wj 1 ' (red line) matches with (WiCj >s Wj r ). This result is obtained from 1000 Monte Carlo simulations 
of our cleaning method using V and W bands. 



2. Cross Power Spectrum 



We now consider point source bias in case of cross power spectrum of two cleaned maps. In this case the cleaned 
power spectrum is given by, 



Cf 



C? + 2(1 - «c)^y+ < WfCJfWf" > + < wfef^w?" > 



(D12) 
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FIG. 17: Left: The point source bias at large multipole range obtained from 1000 Monte Carlo simulations of the cleaning 
method involving V and W bands is shown in blue line. The CMB bias ^7+1 Cf nas b een removed before plotting the blue 
line to show the bias coming from point source only. The red line shows the theoretical power spectrum from which random 
realization of CMB were drawn. Right: Point source bias computed from /(7p !ea ™\ — Cf — 2(1 — n c ) ^-j- is shown in red line. The 

magenta line shows point source bias computed from ^W^Cf^Wf ^> which agrees excellently with W I 1 Cj >s Wf ! computed 
from two arbitrary realizations shown in green and blue lines. 



Here Wj 1 and W 2 are the weight vectors for the two cleaned maps and cf"" 12 '' denotes a chance noise correlation 
in the two sets of maps used to obtain two cleaned maps. From Monte-Carlo simulations involving cross spectrum, 
vlwl2 (g v2w34, vlwl3 ® v2w24, vlwlA <g) v2w23, vlw2A Cg) v2wl3, vlw3A <g) v2w\2 we recover a point source bias at the 
large multipole range. In the left panel of the figure HPTl we show the point source bias from 1000 such simulations. 
Like the auto power case, here too the point source bias could be obtained as W^Cj^Wj 1 . We show this in the right 
panel of the figure [T71 The red line shows bias computed from simulations, i.e. ^Cf lean \ — C\ — 2(1 — n c )i^^. We 
compute Wj^Cf^W 2 for any two arbitrary realizations. These are shown in blue and green colors respectively. The 
magenta colored line shows ^W^C^W 2 \ which is in excellent agreement with blue and green lines. This justifies 

the approximation ^W^C^W 2 \ ~ Wj-C^W 2 , so that point source bias could as well be computed from the 

weights from a single realization without any need of ensemble average. We have followed this approach of point 
source removal in the multi-region iterative cleaning method described in section fill CI 



APPENDIX E: BIAS AT LOW-/ 



Assume that there are nj number of foreground components. Each of these follow rigid frequency scaling. The 
number of channels is n c , where n c > rif + 1. The rank of full covariance matrix is determined by the number of 
independent components. Therefore, for rif number of independent foreground components and an additional CMB 
component the rank of the full covariance matrix is rif + 1. In many cases we choose total number of channels n c to 
be equal to the rank of the matrix, rif + 1. This ensures that the full covariance matrix is of full rank, (ordcr=rank) 
so that it is invertible. 

However it is also interesting to explore the freedom of using a singular covariance matrix Ci. This may arise, for 
example, if we consider n f number of foreground components (each with rigid frequency scaling) but with n c number 
of channels with n c > n/ + 1. The motivation of this work lies in the fact we have shown that the cleaned power 
spectrum could still be defined in terms of Moore Penrose Inversion (in eq. (| 1 7jl ) in case we encounter a singular 
covariance matrix. So it is worth investigating whether there is a further simplified form of eq. (|1T[1 . 

The calculations proceed as follows: 

• Obtain an analytic expression for the full covariance matrix. We find that the full covariance matrix can be 
obtained as three successive rank-one updates of three different matrices. 

• Identify the appropriate cases of the generalized versions of Sherman Morrison formula which applies to each of 
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three rank-one updates that have to be carried out. 

• Apply Sherman Morrison formula successively for the three rank one updates to obtain analytical expression 
for the bias. 

Let the p th foreground component for channel i be denoted by Fq(8, <j>)fp. Here, Fq(9,4>) is the p th foreground 
template based on frequency uo, so that /* = ]., for frequency The full signal map at i th frequency channel is 
given by 



P =i 



Alternatively, in the spherical harmonic space, 



a \m — a lm + E fp' 



pO 

l lm ■ 



(El) 



(E2) 



P =i 



The auto power spectrum of the i th channel 



c\ = df + 2 E f;cfw + E fifpc^ . 



(E3) 



P =i 



p,p> 



where C\ pp '° is the correlation between any two foreground components p,p' and C^^° i s the chance correlation 
between CMB and p th foreground component. The cross power spectrum between two channels i,j are given by 



n } n f n f 

c? = q c + E /pA c/(p)0 + E fi d i mo + E r P f J A pp,) " > 



p=l 



p=l 



or 



p,p' 



c/(p)0 



p— 1 p— 1 

Introducing explicit matrix notations for the equation, we write 

/ 1 1 ... 1 \ / fp fp 



j\ - G; 



1 1 ... 1 



n f 



V 1 • ■•■ 1 / („ c xn 



EA' 
p=i 



■/(p)o 



/p 1 \ 

f 2 f 2 ... f 2 

J p J p J p 



n f 



EA< 
p=i 



■/(p)o 



^ fp fp 
J p J p 



/p c \ 

Jp 



+ 



V/p c 

^ pll pl2 
F 21 p22 



\ fp fp ■■■ fp / 



(n c xn c ) 



\F 



n c l J?n c 2 



J P I (n c xn c ) 
pln c \ 

/ (ll e Xll c ) 



We define 



£pO _ ^c/(p)0 



( ^ \ 
' fl 



V /p" C / 



and 



e = 



1 



(E4) 



(E5) 



(E6) 



(E7) 



(E8) 
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Clearly, f, p0 , eg G R„ c .i, where R nj . )Tle denotes set of real rif x n c matrices. Full covariance matrix can then be written 



as, 



C^QeoeJ+grJeJ + eo^rJ + A3 • (E9) 



Define 



f° = J2 fi P ° . ( E1 °) 
P =i 

then, it is possible to rewrite 

Ci = Q c e eJ + f° e T + e tr + A 3 , (Ell) 

V v ' 

A 2 

V v ' 

Ai 

with f° G M„ c ,i. 



1. Generalized Sherman Morrison Formula for rank one updates 



This section discusses some mathematical theorems that will be useful to us. The results are mainly based upon the 
two papers [H, [29| . In Ref. [28| analytic expressions for MPGI of rank one modified matrices of the form M = A+bc* , 
(notations changed) are reported. Here M, A are any m x n matrices, i.e. M, A G C m ,n, b G C mj i and c G C n 1. 
C mi „ is the set of all m x n complex matrices, (*) denotes conjugate transpose. The motivation of the work in Ref. [28| 
was to generalize the Sherman Morrison formula 

M 1 = A" 1 - ^A- 1 bc*A- 1 (E12) 
A 

where A = 1 + c*A~ 1 b in case of any arbitrary m x n matrix A. The Sherman Morrison formula given by the above 
equation is valid only for square and nonsingular matrix. 

The main results of the Ref. [28[ are a set of formulas for the MPGI depending upon the different set of conditions, 



namely 




(i) 


b £ C(A) and c 


i C(A*) 


(ii) 


b G C(A) and c 


g C(A*) and A = 


(iii) 


b G C(A) and c 


arbitrary and A ^ 


(iv) 


b £ C(A) and c 


G C(A*) and A = 


(v) 


b arbitrary and 


c G C(A*) and A ^ 


(vi) 


b G C(A) and c 


G C(A*) and A = 



where C(A) represents the column space of a matrix A. 

Ref. [2!| shows that it is sufficient to consider 5 independent cases only. All these 5 different cases are listed below 
along with an useful theorem proved regarding the rank modification that occurs during the rank one update of the 
matrix A. 



2. Theorem 



For given A G C m „ and nonzero b G C m> i and c G C-n,i> ^ M be the modifications of A of the form M = A + be* 
and let A = 1 + c*Atb. Then 



r(M) = r(A) - l«be C(A), c g C(A*), A = 
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[ b e C(A),c E C(A*),A ^ 
r(M) = r(A) « be C(A), c ^ C(A*) 
[b £ C(A),c e C(A*) 

r(M) = r(A) + 1 b £ C(A), c ^ C(A*), A = 0. 

Analytical expressions of the MPGI Mt is given in Ref. [lH and in Ref. [29| depending upon the 6 or 5 conditions 
that they find sufficient. The expressions are of the form 

= At + G , (Ef 3) 

where G is a matrix obtained from only sums and products of A, At, b, c and their conjugate transposes. We do not 
mention the explicit form of all the cases here. Instead we give expressions for Mt for those cases which will be useful 
to us. 



a. Case 1 



If b E C(A), c E C(A*), A ^ then 



where 



M 1 = A f - A _1 de* , (E14) 

d = Atb,e= (At)*c. (E15) 

b. Case 2 

If b £ C(A),c i C(A*) then 

M f = At -W - vth + Avt u t, (E16) 
where k = Atb, u = (I AAt)b, v = c*(I - At A), h = c*At. 

c. Case 3 

In Ref. (29[, a relation between the unique projectors between on the column spaces of M, A is reported. If 
b i C(A),c £ C(A*) then 

MM t =AA t -^ 1 ee T + i ) -V 1 qq T , (E17) 
where e = (At)*c, r\ = e*e,v = AA* + rj<p, = f T f , f = (I - AA*)b. 

3. Analytic computation of the bias 

The analytic compution of bias employs the Generalized Sherman-Morrison formula and the relation between 
orthogonal projectors mentioned in the above three cases. Consider the previously stated expression for the covariance 
matrix. 

C, = Q c e eJ + ffej + e fi° T + A 3 . (E18) 



First note that eo ^ C(A3),{° E C(A|), so that this is identical to the conditions for case 3 and rank(As) — 
ranfc(A 2 ) = nt, the number of foreground components. But when we consider matrix A 2 we see that fj° ^ C(A 2 ), e ^ 
C(A 2 ). So here conditions of case 2 apply. Here rank(A 1 ) = rank(A 2 ) + 1 = nt + 1. On the other hand, 
e E C(A!),e E C{A\). Hence rank^^) = ranfc(Ax) = Uf + 1. Now one can carry out the analytic simplification 
in three steps. 



35 



a. Step 1 : First application of GSM 

We note that eo G C(Ai). Therefore we are in a situation where conditions of case 1 are applicable. Our aim is to 
express C^' lean = l/ejcjeo in terms of Ai. For this purpose we simply need to express C[ in terms of Ai following 
GSM formula appropriate for this case. At the next step we compute eJC^eo. The analysis is relatively simple and 
we do not show the detailed calculation here. Instead we give a similar calculation for a more complicated expression 
at the next subsection. The result of this section is, 

Cf lean = Cf + — \r- . (E19) 
e£A[e 



b. Step 2 : Second application of GSM 

At this step we shall simplify l/(ejA^eo) further in terms of some function of A 2 where Ai = f/'ej + A2. As 
mentioned in the previous subsection, this simplification method proceeds in two steps. First, we express A^ in 
terms of A 2 following GSM. At the next stage we would compute l/(ejA] L eo). We easily see that, fj ^ C(A 2 ) and 
eo ^ C(A 2 r ). Therefore GSM formula corresponding to case 2 would be useful for us, 

A{ = A 2 -ku t -v t h + Av t u t . (E20) 



Each of the individual quantities on the right hand side are explained in the section IE 21 All the terms starting from 
the second on the right of above equation could be written in terms of the foreground and CMB shape vectors, fj° 
and eo- The vectors k and h are given by k = A 2 f[° and h = ejA^ respectively. Moore Penrose Inverse of the other 
two vectors u and v could be computed following the definition w = u*/||u||, where ||u|| denotes vector norm. Thus 
we have, 

ut= frd-A 2 At) (I- A tA 2 )eo 
f° T (I-A 2 A 2 )f°' e T(I-A 2 A 2 )e 

All the vectors computed above could be used in eg IE20I to express A^ in terms of A 2 , 

A t = A t _ A 2 f°f° T (I-A 2 A 2 ) _ (I-A 2 A 2 )e e£A 2 T f (I - A 2 A 2 )e f° T (I-A 2 A 2 ) 

2 f° T (I - A 2 At )f,o e T(I-AU 2 )eo ° 2 1 eTfl-AU.leof^d-A^^f, 01 



Now we can easily compute ejA^eo- This involves contraction of both the indices of AJ^.j. After some algebra all 
the terms above simplifies a lot and we are left with, 

1 _ f r(i-A 2 Aj)f» 



ejAle„ f,» T (I - A 2 Ai)e„ 



c. Step 3 : Application of a relation of orthogonal projectors 

Both the numerator and denominator of cq. IE23I contains orthogonal projector A 2 A 2 on the column space of A 2 . 
At this stage we only need to reexpress this term in terms of A3A3 which is also an orthogonal projector on the 
column space of A 3 . We recall that A 2 = A 3 + e fQ. The foreground shape vector fj° lies on the column space 
of the foreground covariance matrix, A3. However the CMB shape vector is linearly independent on the foreground 
templates. Therefore the shape vectors follow, eo ^ C(As), fj° 6 C(A^). Thus the conditions of case 3 are applicable. 
Following the notations of the Ref. [29( we obtain, 

A 2 A 2 = A3A3 - rT 1 ee T + fT 1 ^ 1 qq T (E24) 

To proceed further we need to express ee T and qq T in terms of fj° and eo. We note that, e = A 3 f[° and q = Ae + r/f . 
Following notations of Ref. [29[ we see that f is the component of eo on a plane orthogonal to the column space of 
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A3, i.e. f = (I — A3A 3 )eo. The matrix ee T could immediately be identified as ee T = A 3 fj fj 0T A 3 . However the 
other matrix qq T contains several terms, 

qq T = A 2 ee T + r? 2 (I - A 3 A 3 )e eT(I - A 3 A 3 ) + Xr, (a* ej(l - A 3 A 3 ) + (I - A 3 A 3 )e fp T A 3 ) . (E25) 

It is easy to see that, the numerator of eq. (|E23[) involves computation of the term f 1 0T A2A 2 f 1 °. At this point we 
have all the necessary expressions to compute this term. First we note that, following eqs. (|E24|) . (|E25p the projector 
on the column space of A 2 could be simplified as, 

A2A2 - A3A3 - T^Agfff^Aa + rrV- 1 (A 2 Atf^! 017 A 3 + rf(I - A 3 A 3 )e eJ(I - A 3 A 3 ) + 

Xr, (Atf^e^I - A 3 A 3 ) + (I - A 3 A 3 )e f° T A 3 ) ) . (E26) 



Now we can easily get an expression for f 1 0T A2A 2 f 1 °. Using eq. (|E26|) we obtain, 

f 1 0T A 2 A 2 f 1 ° = f° T A 3 A 3 f° - n-\fr*&?? + V-'v-HXVTKf?) 2 + ^(fi° T (I - A 3 A 3 )e ) 2 + 

2X V (f I 0T A+f 1 °f 1 0T (I - A 3 A 3 )e ) ) . (E27) 

Though there are several terms on the right hand side of the above equation as we will see some of them drop. We 
note that, ff E C(A 3 ), so that (A 3 A 3 )fi° = if. Hence fi° T (I - A 3 A 3 )e = ej(l - A 3 A 3 )fi° = 0. If we assume 
X = f^AsA^fj and Y = f^Agfj then using eq. (|E27|) and v = X 2 + r)(f>, as stated earlier, we may easily obtain, 

f, 0T AaA* f,° = X + ^ . (E28) 



We now identify X = f^AsA^fj and obtain numerator of the eq. (|E23[) as 

f 1 0T (I-A 2 A 2 )f 1 ° = -^. (E29) 



v 



The denominator of eq. (|E23|) can be computed in a similar fashion. We do not elaborate the mathematical details 
any more for this part. The final result is 

f, 0T (I - A 2 A 2 )e = -Y<f>. (E30) 



Thus using eqs. (|E23|) , p28|) , (|E30| we obtain 



ejA[e 

pOT A t f o' 
r i A 3 r l 

simply this term in terms of the rank of the foreground covariance matrix A3. 



' =-Y= -fp T A 3 fi° . (E31) 



The quantity ( T \ — ) — ( f[ Agf" ) constitutes the negative bias in the cleaned power spectrum. Below we further 

e A 1 eo / 



d. Simplification of the bias expression 

We note that the elements of the foreground covariance matrix are given by 

"/ 

^=J2r p f p ,C^°. (E32) 
pp' 

Also the elements of the chance correlation vector between CMB and all the foregrounds 

"/ 

fr=// 0i = EA M °/;- (E33) 
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We may rewrite the magnitude of the negative bias in terms of the above matrix elements and components of the 
foreground shape vector. After a little algebra we get, 

fl 0T At f o^ = £ A W c^C^'A f p f p , . (E34) 

ij \ pp' I 

Using cj- cp )° — J2 m ( a im a, tm)/(^ + -0 an d ^he fact that CMB anisotropies are statistically isotropic, i.e. (a/ m a/' m ') = 
Cf8iv S mm , we may obtain, 

/£ dl^r'A = ^ £ C i PP ' )0 • (E35) 

\ pp' I pp' 

Using eq. (|E35j) and A % £ — J2pp> fpfp'^i PP we can rewrite eq. (|E34j) in terms of following expression consisting of 
CMB and foreground covariance matrices, 

C? 



;f 1 OT Atf 1 o) = -^_^At^. (E36) 
However we observe that, Ylij A$ A$ — tr^A^As) = rank(A3) = rif. Thus we have 

[C? uan ) = (cf)-n } ^-. (E37) 



21 + 1 

which is eq. (|25p. The other formulae given by eq. (j26|) and eq. (|2T[) can be obtained similarly but need much less 
algebra. 
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